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Abstract 

Muons produced in atmospheric cosmic ray showers account for the by far dominant part of the event yield in 
large-volume underground particle detectors. The IceCube detector, with an instrumented volume of about a cubic 
kilometer, has the potential to conduct unique investigations on atmospheric muons by exploiting the large collection 
area and the possibility to track particles over a long distance. Through detailed reconstruction of energy deposition 
along the tracks, the characteristics of muon bundles can be quantified, and individual particles of exceptionally 
high energy identified. The data can then be used to constrain the cosmic ray primary flux and the contribution to 
atmospheric lepton fluxes from prompt decays of short-lived hadrons. 

In this paper, techniques for the extraction of physical measurements from atmospheric muon events are described 
and first results are presented. The multiplicity spectrum of TeV muons in cosmic ray air showers for primaries in 
the energy range from the knee to the ankle is derived and found to be consistent with recent results from surface 
detectors. The single muon energy spectrum is determined up to PeV energies and shows a clear indication for the 
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emergence of a distinct spectral component from prompt decays of short-lived hadrons. The magnitude of the prompt 
flux, which should include a substantial contribution from light vector meson di-muon decays, is consistent with 
current theoretical predictions. 

The variety of measurements and high event statistics can also be exploited for the evaluation of systematic effects. 
In the course of this study, internal inconsistencies in the zenith angle distribution of events were found which indicate 
the presence of an unexplained effect outside the currently applied range of detector systematics. The underlying cause 
could be related to the hadronic interaction models used to describe muon production in air showers. 

Keywords: 

atmospheric muons, cosmic rays, prompt leptons 


1. Introduction 

IceCube is a particle detector with an instrumented 
volume of about one cubic kilometer, located at the ge¬ 
ographic South Pole [T]. The experimental setup con¬ 
sists of 86 cables (“strings”), each supporting 60 digi¬ 
tal optical modules (“DOMs”). Every DOM contains a 
photomultiplier tube and the electronics required to han¬ 
dle data acquisition, digitization and transmission. The 
main active part of the detector is deployed at a depth of 
1450 to 2450 meters below the surface of the ice, which 
in turn lies at an altitude of approximately 2830 meters 
above sea level. The volume detector is supplemented 
by the surface array IceTop, formed by 81 pairs of tanks 
filled with - due to ambient conditions solidified - water. 

The main scientific target of IceCube is the search 
for astrophysical neutrinos. At the time of design, 
the most likely path to discovery was expected to be 
the detection of upward-going tracks caused by Earth- 
penetrating muon neutrinos interacting shortly before 
the detector volume. All DOMs were consequently ori¬ 
ented in the downward direction, such that Cherenkov 
light emission from charged particles along muon tracks 
can be registered after minimal scattering in the sur¬ 
rounding ice. 

The first indication for a neutrino signal exceeding 
the expected background from cosmic ray-induced at¬ 
mospheric fluxes came in the form of two particle show¬ 
ers with a total visible energy of approximately 1 PeV 
El. Detailed analysis of their directionality strongly in¬ 
dicated an origin from above the horizon. The result 
strengthened the case for the astrophysical nature of the 
events, since no accompanying muons were seen, as 
would be expected for neutrinos produced in air show¬ 
ers. This serendipitous detection motivated a dedicated 
search for high-energy neutrinos interacting within the 
detector volume, which led first to a strong indication 
El and later, after evaluating data taken during three 
full years of detector operation, to the first discovery 
of an astrophysical neutrino flux ||4|. In each case, the 


decisive contribution to the event sample were particle 
showers pointing downward. 

Despite the large amount of overhead material, the 
deep IceCube detector is triggered at a rate of approx¬ 
imately 3000 s * by muons produced in cosmic ray- 
induced air showers. Eormerly regarded simply as an 
irksome form of background, these have since proved 
to be an indispensable tool to tag and exclude atmo¬ 
spheric neutrino events in the astrophysical discovery 
region Q. 

Apart from their application in neutrino searches, 
muons can be used for detector verification and a wide 
range of physics analyses. Examples are the measure¬ 
ment of cosmic ray composition and flux in coincidence 
with IceTop la, the first detection of an anisotropy in 
the cosmic ray arrival direction in the southern hemi¬ 
sphere llllISllll, investigation of QCD processes produc¬ 
ing high-pt muons m and the evaluation of track re¬ 
construction accuracy by taking advantage of the shad¬ 
owing of cosmic rays by the moon ini. 

Remaining to be demonstrated is the possibility to de¬ 
velop a comprehensive and consistent picture of atmo¬ 
spheric muon physics in IceCube. The goal of this paper 
is to outline how this could be accomplished, illustrate 
the scientific potential and discuss consequences of the 
actual measurement for the understanding of detector 
systematics. 

2. Physics 

2.7. Cosmic Rays in the IceCube Energy Range 

The energy range of cosmic ray primaries producing 
atmospheric muons in IceCube is limited by the mini¬ 
mum muon energy required to penetrate the ice at the 
low, and the cosmic ray flux rate at the high end. Pre¬ 
dicted event yields are shown in Eig. [T] Since the muon 
energy is related to the energy per nucleon Eprim/A, 
threshold energies increase in proportion to the mass of 
the primary nucleus. 
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Figure 1; Atmospheric muon event yield in IceCube in 
dependence of primary type simulated with CORSIKA 
d. The cosmic ray flux was weighted according to 
the H3a model llT3l . 

The energy range of atmospheric muon events in Ice- 
Cube covers more than six orders of magnitude. Neu¬ 
trinos, not attenuated by the material surrounding the 
detector, can reach even lower. With a ratio between 
lepton and parent nucleon energy of about one order of 
magnitude Gl, the lowest primary energies relevant for 
neutrinos in IceCube fall in the region around 100 GeV. 

Coverage of this vast range of energies by specialized 
detectors varies considerably, and overlapping measure¬ 
ments are not always consistent. At energies well below 
1 TeV, important for production of atmospheric neutri¬ 
nos in oscillation measurements G5l . both PAMELA 
d and AMS-02 Ini And a clear break in the proton 
spectrum at about 200 GeV. The exact behavior of the 
primary spectrum should be an important factor in up¬ 
coming precision measurements of oscillation parame¬ 
ters by the planned IceCube sub-array PINGU ifTSll . 

In the energy region where the bulk of atmospheric 
muons triggering the IceCube detector are produced, 
the most recent measurement was performed by the 
balloon-borne CREAM detector In the range from 
3-200 TeV, proton and helium spectra are found to be 
consistent with power laws of the form The proton 
spectrum with - 2.66 + 0.02 is somewhat softer than 
that of helium with yne = 2.58 + 0.02. The cross-over 
between the two fluxes lies at approximately 10 TeV. 

Between a few hundred GeV and 3 TeV, and again 
from 100 TeV to 1 PeV, there are large gaps where ex¬ 
perimental measurements of individual primary fluxes 
are sparse and contain substantial uncertainties Il20l . Es¬ 
pecially the second region is of high importance to Ice- 
Cube physics, because it corresponds to neutrino ener¬ 
gies of tens of TeV where indications for astrophysical 


fluxes start to become visible. 

The situation improves around the “knee” located at 
about 4 PeV, which has long been a major focus of cos¬ 
mic ray physics. The well constrained overall primary 
flux has been resolved into its individual components by 
the KASCADE array 11211 . although the result depends 
strongly on the model used to describe nuclear interac¬ 
tions within the air shower. There is a general consensus 
that the primary composition changes towards heavier 
elements in the range between the knee and 100 PeV, 
confirmed by various measurements Il22l . including Ice- 
Cube 0. An exact characterization of the all-nucleon 
spectrum around the knee is necessary to constrain the 
contribution to atmospheric lepton fluxes from prompt 
hadron decays and accurately describe backgrounds in 
diffuse astrophysical neutrino searches. 

Between 100 PeV and approximately 1 EeV lies an¬ 
other region with sparse coverage, which has only re¬ 
cently begun to be filled. In the past, data taken near the 
threshold of very large surface arrays indicated a “sec¬ 
ond knee” at about 300 PeV Il2^ . Approaching from 
the other side, KASCADE-Grande found evidence for a 
knee-like structure closer to 100 PeV, along with a hard¬ 
ening of the all-primary spectrum around 15 PeV ll24l . 
This result confirms earlier tentative indications from 
the Tien-Shan detector using data taken before 2001, 
but published only in 2009 ll25l and is supported by 
subsequent measurements using the TUNKA-133 ll26l 
detector. The currently most precise spectrum in terms 
of statistical accuracy and hadronic model dependence 
was derived from data taken by the IceTop surface ar¬ 
ray 1271. KASCADE-Grande later extended the origi¬ 
nal result by indications for a light element ankle ll28l . 
a heavy element knee ll29l and separate spectra for ele¬ 
mental groups lIMIl . 

The emergent picture has yet to be theoretically in¬ 
terpreted in a comprehensive manner. The data indicate 
that several discrete components are present in the cos¬ 
mic ray flux, and that the behavior of individual nuclei 
closely corresponds to a power law followed by a spec¬ 
tral cutoff at an energy proportional to their magnetic 
rigidity R - This explanation was first pro¬ 

posed by Peters in 1961 lITTI and later elaborated by, 
among others, Ter-Antonyan and Haroyan 021 as well 
as Horandel 03l . Exactly how many components there 
are, where they originate, and the precise values and 
functional dependence of their transition energies are 
still open questions. A well-known proposal by Hillas 
postulates two galactic sources, one accounting for the 
knee, the other for the presumptive knee-like feature at 
300 PeV 0^ . Another model, by Zatsepin and Sokol- 
skaya, identifies three distinct types of galactic sources 
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to account for the flux up to 100 PeV lll5]l . The harden¬ 
ing of the spectrum around the “ankle” at several EeV 
can be described elegantly by a pure protonic flux and 
its interaction with CMB radiation 1361 or, more in line 
with recent experimental results, in terms of separate 
light and heavy components l37l . The consensus is in 
either case that the origin of the highest-energy cosmic 
rays is extragalactic. 

This paper, like other IceCube analyses, relies for 
purposes of model testing mainly on the parametriza- 
tions by Gaisser, Stanev and Tilav na. These incor¬ 
porate various basic features of the models described 
above, while updating numerical values to conform with 
the latest available measurements. Specifically, the 
“Global Fit” (GF) parametrization introduces a second 
distinct population of cosmic rays before the knee with 
a transition energy of 120 TeV. The knee itself, and the 
feature at 100 PeV, are interpreted as helium and iron 
components with a common rigidity-dependent cutoff, 
eliminating the need for an intermediate galactic flux 
component as in the H(illas) 3 a and 4a parametrizations. 
The difference between H3a and H4a lies in the compo¬ 
sition of the highest-energy component which becomes 
dominant at energies beyond 1 FeV, which is mixed in 
the former, and purely protonic in the latter case. In the 
region around the knee, the two models are for practical 
purposes indistinguishable. 

2.2. Muons vs. Neutrinos 

The flux of atmospheric neutrinos in IceCube is mod¬ 
eled using extrapolated parametrizations based on a 
Monte Carlo simulation for energies up to 10 TeV ||39|. 
To account for the influence of uncertainties of the cos¬ 
mic ray nucleon flux, the energy spectrum is adjusted by 
a correction factor BOl . The result can be demonstrated 
to agree reasonably well with full air shower simula¬ 
tions m, but necessarily contains inaccuracies, for ex¬ 
ample by neglecting variations in the atmospheric den¬ 
sity profile at the site and time of production. 

Atmospheric muon events on the other hand are sim¬ 
ulated through detailed modeling of individual cosmic 
ray-induced air showers. In standard simulation pack¬ 
ages such as CORSIKA QIl, specific local conditions 
like the direction of the magnetic field and the profile 
of the atmosphere including seasonal variations can be 
fully taken into account. Energy spectra for each type of 
primary nucleus are separately adjustable. Hadronic in¬ 
teraction models can be varied and their influence quan¬ 
tified in terms of a systematic uncertainty. 

In contrast to neutrinos, astrophysical fluxes, flavor¬ 
changing effects and hypothetical exotic phenomena do 


not affect muons. All observations can be directly re¬ 
lated to the primary cosmic ray flux and the detailed 
mechanisms of hadron collisions. Due to the close re¬ 
lation between neutrino and charged lepton production, 
high-statistics measurements using muon data are there¬ 
fore invaluable to constrain atmospheric neutrino fluxes. 

Perhaps most importantly, atmospheric muons repre¬ 
sent a high-quality test beam for the verification of de¬ 
tector performance, because the variety of possible mea¬ 
surements along with high event statistics permit de¬ 
tailed consistency checks. A particular advantage in the 
case of IceCube is that muons probe the region above 
the horizon for which the down-looking detector con¬ 
figuration is not ideal, but where contrary to original ex¬ 
pectation the bulk of astrophysical detections has taken 
place. 

2.3. Primary Flux and Atmospheric Muon Character¬ 
istics 

The connection between the measurable quantities of 
atmospheric muon events and the properties of the pri¬ 
mary cosmic ray flux is illustrated in Fig. The rela¬ 
tion of muon multiplicity to primary type and energy 
is expressed in terms of the parameter Fmuit, defined 
such that Fmuit = Fprim for iron primaries. The average 
number of muons in a bundle can then be expressed as 
< N/t >- K ■ Frnuit, where the proportionality factor k de¬ 
pends on the specific experimental circumstances. Due 
to fluctuations in the atmospheric depth of shower de¬ 
velopment and the total amount of hadrons produced in 
nuclear collisions, the variation in the number of muons 
is slightly wider than a Poissonian distribution IMll . 

Since the muon multiplicity itself is a function of 
zenith angle, atmospheric conditions, detector depth 
and surrounding material, it is convenient to re-scale it 
such that the derived quantity is directly related to pri¬ 
mary mass and energy. This study uses the parameter 

^mult = -Epi-im ■ (A/56)^. (1) 

The definition was chosen such that Fmuit is equal to 
Fprim in the case of iron primaries with atomic mass 
number A - 56, which will in practice dominate the 
multiplicity spectrum above a few PeV, as shown in Fig. 
[^(b). Exact definition and construction of Fmuit are dis¬ 
cussed in Section lhTI 

As the ratio of muons to electromagnetic particles in 
an air shower increases with the primary mass, the con¬ 
tribution of light elements to the multiplicity spectrum is 
suppressed. For a power law spectrum of the form E^'^, 
the contribution of individual elements to the muon mul¬ 
tiplicity, here expressed in terms of a flux Chmuit, scales 
as: 
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(a) Primary Energy 


where a « 0.79 is an empirical parameter derived from 
simulation m. 

Single-particle atmospheric lepton fluxes, on the 
other hand, are related to the nucleon spectrum. Un¬ 
der the same assumptions as above, the relation between 
all-nucleon and primary flux is: 


‘I*nuc(^) *^prim(l) _ 

<l>nuc(l) ‘tprimCA) 

For a power law with an index of approximately -2.6 to 
-2.7, such as the cosmic ray spectrum before the knee, 
the nucleon spectrum therefore becomes strongly dom¬ 
inated by light elements. 



(b) Rescaled Muon Multiplicity Enmn 



(c) Nucleon Energy 


Figure 2: Contribution of individual elemental compo¬ 
nents to overall flux spectra relevant for atmospheric 
muon measurements, here shown for the Gaisser/Hillas 
model with mixed-composition extragalactic compo¬ 
nent (H3a) 113. For definition of Tsmuit, see Section|2.3 


^prim(l) 

^prim(.'^) 


2.4. Prompt Muon Production 

A particular difficulty in the description of atmo¬ 
spheric lepton fluxes is the emergence at high energies 
of a component originating from prompt hadron decays. 
The reason is the harder spectrum compared to the light 
meson contribution, which is the consequence of the 
lack of re-interactions implicit in the definition. 

An important source of prompt atmospheric lepton 
fluxes is the decay of charmed hadrons. While it is 
possible to estimate their production cross section us¬ 
ing theoretical calculations based on perturbative QCD, 
substantial contributions from non-perturbative mecha¬ 
nisms cannot be excluded. The problem can therefore at 
the moment only be resolved experimentally B2ll . One 
major open question currently under investigation 143 
is whether nucleons contain “Intrinsic Charm” quarks, 
which might considerably increase charmed hadron pro¬ 
duction B4ll . 

Inclusive charm production cross sections were mea¬ 
sured during recent LFIC runs by the collider detectors 
LHCb 113, ATLAS 113, and ALICE l47ll48l . and pre¬ 
viously by the RHIC collaborations PHENIX l49l and 
STAR l50l . Data points are consistently located at the 
upper end of the theoretical uncertainty, which covers 
about an order of magnitude ISTIl . On a qualitative level, 
the new results suggest that charm-induced atmospheric 
neutrino fluxes could be somewhat stronger than previ¬ 
ously assumed. A straightforward translation is how¬ 
ever difficult. Although collider measurements probe 
similar center-of-mass energies, they are for technical 
reasons restricted to central rapidities of approximately 
lyl < 1. Eor lepton production in cosmic ray interac¬ 
tions, forward production is much more important. 

A variety of descriptions for the flux of atmospheric 
leptons from charm have been proposed in the past 153 . 
In recent years, the model by Enberg, Reno and Sarce- 
vic 153 has become the standard, especially within the 
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IceCube collaboration, which usually expresses prompt 
fluxes in “ERS units”. For muons, electromagnetic de¬ 
cays of unflavored vector mesons make a significant ad¬ 
ditional contribution not present in neutrinos ll54l , and 
at the very highest energies di-muon pairs are produced 
by Drell-Yan processes Il55l . Especially the first process 
should lead to a substantial enhancement of the prompt 
muon flux compared to neutrinos ll5^ . A detailed dis¬ 
cussion can be found in Appendix B| 

It has long been suggested to use large-volume neu¬ 
trino detectors to constrain the prompt component of the 
atmospheric muon flux directly lISTl . Apart from the 
aspect of particle physics, the approximate equivalence 
between prompt muon and neutrino fluxes would help to 
constrain atmospheric background in the energy region 
critical for astrophysical searches. 

Past measurements of the muon energy spectrum in 
volume detectors were not able to identify the prompt 
component. Usually based on the zenith angle distribu¬ 
tion alone, the upper end of their energy range fell one 
order of magnitude or more below the region where the 
prompt flux is expected to become measurable ll58l . The 
LVD collaboration, by exploiting azimuthal variations 
in the density of the surrounding material, was able to 
set a weak limit ll59l . The Baksan Underground Scin¬ 
tillation Telescope reported a significant excess above 
even the most optimistic predictions fbO), but the result 
has not yet been confirmed independently. 


3. Data Samples 

3.1. Experimental Data 

The data used in this study were taken during two 
years of detector operation from 2010 to 2012. Origi¬ 
nally the analysis was developed for the first year only, 
but problems related to simulation methods as discussed 
in Section [T2| made it necessary to base the high-energy 
muon measurement on the subsequent year instead. 


Time Period 

Detector Configuration 

Livetime 

05-31-2010-05-13-2011 
05-13-2011 -05-15-2012 

79 Strings (IC79) 

86 Strings (IC86) 

313.3 days 
332.1 days 


Table 1: Experimental Data Sets. 

The main IceCube trigger requires four or more pairs 
of neighboring or next-to-neighboring DOMs to regis¬ 
ter a signal within a time of 5 ps. Full event information 
is read out for a window extending from 10 ps before 
to 22 ps after the moment at which the condition was 
fulfilled. Including events triggered by the surface ar¬ 
ray IceTop and the low-energy extension DeepCore, for 


which special conditions are implemented, this results 
in a total event rate of approximately 3000 s ' for the 
full 86-string detector configuration. 

As data transfer from the South Pole is constrained by 
bandwidth limitations, only specific subsets are avail¬ 
able for offline analyses. The main requirement in the 
studies presented here was an unbiased base sample. 
The main physics analyses therefore use the filter stream 
containing all events with a total of more than 1,000 
photo-electrons. Additionaly, minimum bias data corre¬ 
sponding to every 600th trigger were applied to evaluate 
detector systematics. 

Reconstruction of track direction and quality parame¬ 
ters followed the standard IceCube procedure for muon 
candidate events 16^ . based on multiple photo-electron 
information and including isolated DOMs registering a 
signal. In addition, various specific energy reconstruc¬ 
tion algorithms were applied. For all data, the differ¬ 
ential energy deposition was calculated using the de¬ 
terministic method discussed in Appendix A and the 


track energy was estimated by a truncation method | 
Likelihood-based energy reconstructions lf64l were ap¬ 
plied to the first year of data only, primarily for evalua¬ 
tion purposes. 


3.2. Simulation 

The standard method used for simulation of cosmic 
ray-induced air showers in IceCube is the CORSIKA 
software package IfT^ . in which the physics of hadronic 
interactions are implemented via externally developed 
and freely interchangeable modules. In this study, as in 
all IceCube analyses, mass air shower simulation pro¬ 
duction was based on SIBYLL 2.1 ll65ll . To investigate 
systematic variations, smaller sets of simulated data 
were produced using the QGSJET-II and EPOS 
1.99 El models. 

In the cuiTent version of CORSIKA (7.4), the con¬ 
tribution from prompt decays of charmed hadrons and 
short-lived vector mesons to the muon flux is usually 
neglected. An accurate simulation would in any case 
be difficult due to strong uncertainies on production 
and re-interaction cross sections. For this study, the 
prompt component of the atmospheric muon flux was 
modeled by re-weighting events produced in decays of 


light mesons. The exact procedure is described in Ap- 
pendix B| 


High-energy muons passing through matter lose their 
energy through a variety of specific processes ll68l . 
which in IceCube are modeled by a dedicated simula¬ 
tion code l69l . The energy spectra of discrete catas¬ 
trophic losses along atmospheric muon tracks predicted 
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Figure 3: Energy spectra of discrete stochastic energy 
losses along muon tracks simulated using the mmc code 
1^ . The data sample corresponds to events with more 
than 1,000 registered photo-electrons in the IceCube de¬ 
tector. For demonstration purposes, the primary cosmic 
ray spectrum is modeled as an unbroken ^ power 
law. 

to occur within the IceCube detector volume are shown 
in Fig. 

For all energy loss processes, the corresponding 
Cherenkov photon emission is calculated. Every pho¬ 
ton is then tracked through the detector medium until it 
is either lost due to absorption or intersects with an op¬ 
tical module llTOl . This detailed procedure is necessary 
to account for geometrically complex variations in the 
optical properties of the ice, but has the disadvantage of 
being computationally intensive, limiting the amount of 
simulated data especially for bright events. 

The variations between direct photon propagation 
and the tabular method previously used in IceCube sim¬ 
ulations were evaluated for each of the studies presented 
in this paper. It was found that in the case of high- 
multiplicity bundles the difference can be accounted for 
by a simple correction factor, while for high-energy 
tracks the distortion was so severe that simulations pro¬ 
duced with the obsolete method were unusable. Simu¬ 
lation mass production based on direct photon propaga¬ 
tion is only available for the 86-string detector configu¬ 
ration, requiring the use of a corresponding experimen¬ 
tal data set. In order to reduce computational require¬ 
ments, the measurement of bundle multiplicity was not 
duplicated and instead solely relies on data from the 79- 
string configuration. 

The low cosmic ray flux rate at the highest primary 
energies means that even relatively few events corre¬ 
spond to large amounts of equivalent live time. Accord¬ 
ingly, for the measurement of the bundle multiplicity 


spectrum simulation statistics are not a limiting factor. 
In the region before and at the knee, where the domi¬ 
nant part of high-energy muons are produced, far more 
showers need to be simulated. For this reason, the statis¬ 
tical accuracy of the single muon energy spectrum mea¬ 
surement is limited by the amount of simulated livetime, 
generally corresponding to substantially less than one 
year. 

The calculation of detector acceptance and conver¬ 
sion of muon fluxes from South Pole to standard condi¬ 
tions for high energy muons as described in Section [74| 
made use of an external simulated data set produced for 
a dedicated study on the effect of hadronic interaction 
models on atmospheric lepton fluxes ED. 

4. Low-Energy Muons 

4.1. Observables 

A comprehensive verification of detector perfor¬ 
mance requires the demonstration that atmospheric 
muon data are understood at a basic level. Sufficient 
statistics for this purpose are in IceCube provided by the 
minimum bias sample, consisting of every 600th event 
triggering the detector. 

Two simple parameters were used in the evaluation. 
These are the zenith angle of the reconstructed 
track, with flzen = 0 for vertically down-going muons 
from zenith, and the total number of photo-electrons 
2tot registered in the event. 

The angular dependence of the muon flux can be di¬ 
rectly related to the energy spectrum in the TeV range, 
because the threshold increases as a function of the 
amount of matter that a muon has to traverse before 
reaching the detector. The limiting factors near the hori¬ 
zon are the rapid increase of the mean surface energy 
approximately proportional to exp(sec 0zen)j the corre¬ 
sponding decrease in flux, and eventually the irreducible 
background from atmospheric muon neutrinos. Purely 
angular-based muon energy spectra therefore only reach 
up to energies of 20-30 TeV, depending on the depth of 
the detector and the type of surrounding material. For 
the specific case of IceCube, the relation of zenith angle 
to muon and primary nucleon energy is shown in Fig. 

The total number of photo-electrons (“brightness”) 
of atmospheric muon events is closely related to the 
muon multiplicity, as demonstrated in Fig. where 
events with photons registered by the DeepCore array 
were excluded to avoid minor biases at the very low 
end of the distribution. In the experimental measure¬ 
ments described below, all events were included. The 
emitted Cherenkov light is in good approximation pro¬ 
portional to the total energy loss, and the multiplicity 
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(a) Muon Surface Energy 



log 


10 "’ 


(b) Parent Nucleon Energy 


Figure 4; Relation between reconstructed zenith angle 
and energy for simulated muon showers triggering the 
IceCube detector. The distributions correspond to min¬ 
imum bias data after track quality selection described 
in Sec. 4.3 Superimposed are mean and spread of the 
distribution. 



log,„Q,„/[pe] 


Figure 5; Top; Simulated distributions of total number 
of photo-electrons in event, separated in dependence of 
number of muons in bundle at closest approach to the 
center of the IceCube detector. The functional depen¬ 
dence of the fit is described in the text. Bottom: Change 
of data/simulation ratio for different assumptions about 
the light yield, effectively corresponding to the relation 
between energy deposition and number of registered 
photo-electrons. The simulation was weighted accord¬ 
ing to the H3a primary flux model. 


spectrum can therefore be measured even at low ener¬ 
gies, although its interpretation is difficult because of 
the varying threshold for the individual components of 
the cosmic ray flux. 

The distribution for a hxed number of muons can 
be described by a transition from a Gaussian distribu¬ 
tion to an exponential in terms of the parameter q = 
logio(2tot/p-e.); 

^Wevent ,, ~ ~ Pfi{q — q^eak) 

-:- =A'-eXp -^ H--r 

I 1 + ^peak) \ + exp ?pcak) 

(4) 

The free ht parameters for the case of single muon 
events are described in Table |2] While all values de¬ 
pend on the exact detector setup and event sample and 


Fit Parameter 

Value 

Interpretation 

^peak 

1.615 ±0.002 

42.2 p.e 

a 

5.35 ±0.34 

Transition Smoothness 

(T 

0.160 ±0.004 

Width of Gaussian 

P, 

-6.23 ± 0.07 

Power Law Index 

N 

arbitrary 

normalization 


Table 2: Parameters and values for the fit to the single 
muon distribution shown in Fig. i The ;if^/dof of the 
ht is 26.75/16, where the main deviation from the ht is 
found in the hrst three bins of the histogram. 


have no profound physical meaning, the description 
nevertheless provides valuable insights. For example. 
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the peak position corresponds to the average number 
of photo-electrons detected from a minimum ionizing 
track crossing the full length of the detector, and repre¬ 
sents an approximate calorimetric scale from which the 
response to a given energy deposition can be estimated. 

The lower, Gaussian half of the one-muon distribu¬ 
tion only depends on the experimental setup and shows 
minimal sensitivity to physics effects in simulations. In 
particular, the peak value g'peak varies as a function of the 
optical efficiency, a scalar parameter which expresses 
the effects of a wide variety of underlying phenomena 
EH. As shown in the lower panel of Fig. above a 
certain threshold only the flux level, not the shape of 
the distribution is affected by detector systematics. This 
is a common observation for energy-related observables 
and a simple consequence of the effect of a slight offset 
on a power law function. Note that the measured dis¬ 
tribution is fully consistent with expectation within the 
10% light yield variation usually assumed as systematic 
uncertainty in IceCube. 

4.2. Connection to Primary Flux 

The consistency of measurements on separate observ¬ 
ables can be checked by relating them to the primary 
cosmic ray flux. Assuming that the current understand¬ 
ing of muon production in air showers is correct, there 
should be a model which describes both energy and 
multiplicity spectra of atmospheric muons. 

Figure shows the two proxy variables described in 
the previous section, separated by elemental type of the 
cosmic ray primary. At all angles, the muon flux is 
strongly dominated by proton primaries. This is a sim¬ 
ple consequence of the connection between muon en¬ 
ergy and energy per nucleon of the primary particle, and 
does not depend strongly on the specific cosmic ray flux 
model ED- Likewise, the multiplicity-related bright¬ 
ness distribution is for low values dominated by light 
primaries, a consequence of the varying threshold ener¬ 
gies shown in Fig. [T] 

The cosmic ray flux models best reproducing the lat¬ 
est direct measurement in the relevant energy region 
from 10 to 100 TeV HD are GST-GF and H3a US. For 
the comparisons between data and simulation in the fol¬ 
lowing section, they are used as benchmark models rep¬ 
resenting the best prediction at the current time. In ad¬ 
dition, toy models corresponding to straight power law 
spectra are discussed to illustrate the effect of variations 
in the primary nucleon index. In these, elemental com¬ 
position and absolute flux levels at 10 TeV primary en¬ 
ergy correspond to the rigidity-dependent poly-gonato 
model ll^ . used as default setting for the production of 
IceCube atmospheric muon systematics data sets. 




Figure 6: Low-level observables for IceCube atmo¬ 
spheric muon events at trigger level, separated by cos¬ 
mic ray primary type. The simulated data were gener¬ 
ated with CORSIKA IfT^ and weighted according to the 
H3a model IIT3l . 

4.3. Experimental Result 

For the study presented in this section, minimum bias 
data and simulation were compared at trigger level and 
for a sample of high-quality tracks requiring: 

• Reconstructed track length within the detector ex¬ 
ceeding 600 meters. 

• llhiecoKNuoM - 2.5) < 7.5, where llh^eco corre¬ 
sponds to the likelihood value of the track recon¬ 
struction and Vdom to the number of optical mod¬ 
ules registering a signal. 

The stringency of the quality selection is slightly 
weaker than in typical neutrino analyses. For tracks re¬ 
constructed as originating from below the horizon, the 
contribution from mis-reconstructed atmospheric muon 
events amounts to about 50%. 

Simulated and experimental zenith angle distribu¬ 
tions are shown in Fig. |7] Even at trigger level, the influ¬ 
ence of mis-reconstructed tracks can be neglected in the 
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Figure 7: Angular distribution of true and reconstructed 
atmospheric muon tracks in simulation compared to ex¬ 
perimental data. Top: Trigger Level, Bottom: High- 
Quality Selection. The event sample corresponds to 
minimum bias data encompassing all trigger types. The 
ratio of experimental data to simulation is shown in 
Figs. |^(a) and (c). 


region above 30 degrees from the horizon (cos flzen = 
0.5). For the high-quality data set, true and recon¬ 
structed distributions are approximately equal down to 
angles of cos 0zen = 0.15, or 80 degrees from zenith. 

Figure]^ shows comparisons between data and simu¬ 
lation weighted according to several primary flux pre¬ 
dictions. The total number of photo-electrons is de¬ 
scribed reasonably well by the simulation weighted ac¬ 
cording to the H3a model. Application of quality cri¬ 
teria does not lead to any visible distortion. The an¬ 
gular distribution, on the other hand, shows substantial 
inconsistencies. At trigger level, the spectrum is clearly 
harder than for the high-quality sample. The discrep¬ 
ancy does not depend on the particular track quality pa¬ 
rameters used in the selection. 

It is important to note that the absolute level of the 
ratio is not a relevant quantity for the evaluation. Con¬ 


sistency between measurement and expectation within 
the range of systematic uncertainties on the photon yield 
was demonstrated for the brightness distribution in Sec¬ 
tion 4.1 Also, absolute primary flux levels derived from 


direct measurements are typically constrained no better 
than to several tens of percents. For the toy models, 
the normalization was consciously chosen to produce a 
clear separation from the realistic curves in the interest 
of clarity. 

The trigger-level angular distribution in the re¬ 
gion near the horizon becomes dominated by mis- 
reconstructed events consisting of two separate showers 
crossing the detector in close succession. The frequency 
of these “coincident” events scales with the square of 
the overall shower rate, leading to a spurious distortion 
of the ratio between data and simulation in cases where 
the absolute normalization is not exactly equal. This ef¬ 
fect is visible in Fig. [^(a) at values below 0.3. 

To quantify the discrepancy between trigger and 
high-quality level and investigate the influence of sys¬ 
tematic uncertainties, the toy model simulation was fit¬ 
ted to data for 1 > cos flzen > 0.5. In this region, in¬ 
fluences of mis-reconstructed tracks are negligible even 
at trigger level, as demonstrated in Fig. |7] From Fig. 
it can be seen that this corresponds to a relatively small 
energy range for muons and parent nuclei, over which 
the power law index of the cosmic ray all-nucleon spec¬ 
trum can be assumed to be approximately constant and 
used as sole fit parameter. As the normalization was left 
free, the best result simply corresponds to a flat curve 
for the ratio between data and simulation. Possible ef¬ 
fects of variations in the primary elemental composition 
can be taken into account as a systematic error. 

The numerical results of the fit to the angular distribu¬ 
tion is shown in Tabled Note that for cases where the 
statistical error due to limited simulated data exceeds 
the absolute value of the variation, only an upper limit is 
given. The best fit results for the spectral index at trigger 
and high-quality level, 2.715 and 2.855, are illustrated 
by the toy model curves in Fig. Both measurements 
are softer than those of the realistic models, in which 
Tnucleon ~ 2.64. 


4.4. Interpretation 

For the strong discrepancy between the measure¬ 
ments at trigger and high-quality level of AycR = 
0.140 + 0.008(stat.), the following explanations can be 
proposed: 

• A global adjustment to the bulk ice absorption 
length of more than 20%. This explanation would 
imply a major flaw in the method used to derive the 
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(a) cos Szen. Trigger Level 




(c) cos 0zein High-Quality Tracks 



(d) Qtot, High-Quality Tracks 


Figure 8; Ratio of experimental data to simulation in terms of reconstructed zenith angle 0zen and total amount of 
registered photo-electrons Qtot- The primary flux models used in this comparison are discussed in Section 


4.2 


Type 

Variation 

yCR.Tiigger 

yCR,High-Q 

AycR 

Hole Ice Scattering 
Bulk Ice Absorption 
Bulk Ice Scattering 
Primary Composition 
Hadronic Model 
DOM Efficiency 
Experimental Value 

30cm/100cm 
+ 10% 

+ 10% 
p/He 

QGSJET-II/EPOS 1.99 
+ 10% 

Statistical Error 

+0.03 

+0.03 

<0.01 

<0.01 

+0.02/ < 0.01 
<0.02 

2.715 + 0.003 

+0.03/ - 0.05 
+0.02 
+0.01 

+0.03/-0.10 
+0.03/ < 0.02 
+ < 0.02/ - 0.04 
2.855 + 0.007 

+0.01/-0.02 
+0.05 
<0.015 
-0.03/+ 0.10 
< 0.02 

+0.02/- < 0.02 
0.140 + 0.008 


Table 3; Cosmic ray nucleon spectrum measurement and influence of systematic uncertainties. The goodness of the 
experimental fit is ^If^/dof = 13.0/11 at trigger and 12.6/11 at high-quality level. 


optical ice properties Il72l . and is strongly disfa¬ 
vored by the good agreement between the effective 
attenuation length in data and simulation demon¬ 


strated in Appendix A 


• A substantial change of the primary cosmic ray 
composition towards heavier elements. In an event 


sample entirely excluding proton primaries, the 
observed effect can be approximately reproduced. 
However, the increased threshold energy would re¬ 
quire the overall primary flux to be more than three 
times higher than in the default assumption to pro¬ 
duce the observed event rate. An explanation based 
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purely on a heavier cosmic ray composition there¬ 
fore appears highly unlikely. 

• A major inaccuracy of hadronic interaction simula¬ 
tions common to SIBYLL, QGSJET-II and EPOS. 
While this explanation seems improbable, espe¬ 
cially given the almost perfect agreement between 
SIBYLL and EPOS, it should be noted that the 
models used in the IceCube CORSIKA simulation 
were developed before LHC data became avail¬ 
able. Improved models are in preparation ll7^ and 
it should be possible to evaluate them in the near 
future. 

• An unsimulated detector effect with a significant 
influence on the behavior of track quality parame¬ 
ters. Detectors using naturally grown ice are inher¬ 
ently difficult to model in simulations. The optical 
properties of the medium are inhomogeneous and 
photon scattering has a substantial influence on the 
data. The situation is complicated further by the 
placement of the active elements in re-frozen “hole 
ice” columns containg sizable amounts of air bub¬ 
bles. Studies on possible error sources are ongoing 
at the time of writing, but currently there is no in¬ 
dication for a major oversight. 

While the presence of an inconsistency is clear, from 
IceCube data alone there is no strict way to conclude 
whether the brightness or the angular measurement is 
more reliable. However, the evidence strongly points 
to an unrecognized angular-dependent effect introduced 
by track quality-related observables. The reasons are: 

• The brightness distributions are consistent both be¬ 
tween the two data samples and with direct mea¬ 
surements of the cosmic ray flux. 

• At trigger level, there are no major contradictions 
between brightness and zenith angle distributions. 

• The angular spectrum for the high-quality data 
set is significantly steeper than both the neutrino- 
derived result lED and direct measurements. In 
comparisons to the latter, the error from the varia¬ 
tion in primary composition does not apply, as pro¬ 
ton and helium fluxes are constrained individually. 
The total systematic uncertainty on the all-nucleon 
power law index would in this case be reduced to 
about +0.06, whereas the difference in measure¬ 
ment is larger than 0.2. On the other hand, it is 
interesting to note that the LVD detector found a 
value of jcr - 2.78 + 0.05 ll74l . closer to the Ice- 
Cube high-quality sample result. 




Eigure 9: Distribution of muon energies in individual air 
showers at the IceCube detector depth simulated with 
CORSIKA/SIBYLL 112ISl, averaged over all angles. 
Top: Eprim = 3 PeV. Bottom: Epiim = 100 PeV. The 
threshold effect visible at high muon energies in the 
top plot is due to the lower energy per nucleon in iron 
primaries. As the total energy increases, this effect be¬ 
comes less and less visible and the spectra are identical 
except for a scaling factor. 


Even though angular distributions of atmospheric 
muons have been published by practically all large- 
volume neutrino detectors and prototypes ll75l 17^ l77l 
|28]|29l[80]|8Tl[82l, none of the measurements is accu¬ 
rate enough to provide a strict external constraint. Eor 
the time being, there is no other choice than to note the 
effect and continue to investigate possible explanations. 
In the main physics analyses described in the subse¬ 
quent sections, the possible presence of an angular dis¬ 
tortion was taken into account as a systematic error on 
the result. 
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5. Physics Analyses 


While the study of low-energy atmospheric muons is 
instructive for detector verification and the evaluation of 
systematic uncertainties, the main physics potential lies 
in the measurement of events at higher energies. Here it 
is necessary to distinguish two main categories: 

• High-Multiplicity Bundles, in which muons con¬ 
form to typical energy distributions as shown in 
Fig. 0 The total energy 2 contained in the 
bundle is approximately proportional to the num¬ 
ber of muons and related to primary mass A 
and energy Fprim as 


(5) 

with a =» 0.79. The dependence of the muon mul¬ 
tiplicity on the mass of the cosmic ray primary is 
the main principle underlying composition analy¬ 
ses using deep detector and surface array in coin¬ 
cidence 161. Low-energy muons lose their energy 
smoothly, and fluctuations in the energy deposition 
are usually negligible. 

• High-Energy Muons with energies significantly 
exceeding the main bundle distribution. Their pro¬ 
duction is dominated by exceptionally quick de¬ 
cays of pions and kaons at an early stage in the de¬ 
velopment of the air shower. Figure shows that 
showers with more than one muon with an energy 
above several tens of TeV are very rare. Any muon 
with an energy of 30 TeV or more will therefore 
very likely be the leading one in the shower, al¬ 
though this does not exclude the presence of other 
muons at lower energies. The primary nucleus can 
in this case be approximated as a superposition of 
individual nucleons, each carrying an energy of 
Enucieon = Eprim/A. High-energy lepton spectra are 
therefore a function of the primary nucleon flux. 

Hadronic models, cosmic ray spectrum and com¬ 
position all have a significant influence on TeV 
muons ll83i . In addition, at muon energies ap¬ 
proaching 1 PeV prompt decays of short-lived 
hadrons play a significant role. The result is a com¬ 
plex picture with substantial uncertainties, as nei¬ 
ther the exact behavior of the nucleon spectrum at 
the knee nor the production of heavy quarks in air 
showers is fully understood. A schematic illustra¬ 
tion of the muon flux above 100 TeV is given in 
Fig.fTT 



Figure 10: Surface energy distribution for all and most 
energetic (“leading”) muons in simulated events with a 
total of more than 1,000 registered photo-electrons in 
IceCube. 


Charged leptons and neutrinos are usually pro¬ 
duced in the same hadron decay. The energy spec¬ 
trum of single muons is therefore the quantity most 
relevant for the constraint of atmospheric neutrino 
fluxes. Since the stochasticity of energy losses in 
matter increases with the muon energy, the signal 
registered in the detector can vary substantially, as 
in the case of neutrino-induced muons. 



Figure 11: Sketch illustrating the contribution to the sin¬ 
gle muon spectrum at energies beyond 100 TeV. The 
“conventional” component from light mesons is sensi¬ 
tive to atmospheric density and varies as a function of 
the zenith angle JH, that from prompt decays of short¬ 
lived hadrons is isotropic. Re-interactions cause the 
non-prompt spectrum to be steeper. The exact spectral 
shape depends on the all-nucleon cosmic ray flux, with 
a significant steepening expected due to the cutoff at the 
“knee”. 
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The transition between the two atmospheric muon 
event types is gradual. High-energy events rarely con¬ 
sist of single particles, and the characteristics of the ac¬ 
companying bundle of low-energy muons could in prin¬ 
ciple for some cases be determined and used to extract 
additional information about the primary nucleus. At 
low energies the distinction becomes meaningless, as 
events are usually caused by single or very few muons 
with energies below 1 TeV. 
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Figure 12: Event samples used for the measurements 
described in Sec. |^and|7] Shown are true parameter 
distributions for simulated data with more than 1,000 
registered photo-electrons. Top: Fraction of total bundle 
energy carried by the leading muon. Bottom: Energy 
of CR primary. The bimodal shape of the distributions 
becomes more pronounced with increasing brightness. 


Two separate analysis samples were extracted from 
the data, corresponding to high-energy muon and high- 
multiplcity bundle event types. Figure 12 illustrates 
their characteristics in terms of true event parameters 
derived from Monte-Carlo simulations. High-energy 
events, in which the total muon energy is dominated by 
the leading particle, are outnumbered by a factor of ap¬ 
proximately ten. The corresponding need for more rig¬ 


orous background suppression leads to a lower selection 
efficiency than in the case of large bundles. The details 
of the selection methods are described in the following 
sections. 

6. Muon Bundle Multiplicity Spectrum 

6.1. Principle 

The altitude of air shower development, and with 
it the fraction of primary energy going to muons, de¬ 
creases as a function of parent energy Eprim, but in¬ 
creases with the nuclear mass A. The relation between 
the energy of the cosmic ray primary and the number 
of muons above a given energy is therefore not 

linear. A good approximation is given by the “Elbert 
formula”: 


Ap(£ > = A- 


P-ii.mm. COS d 


AE„ 


1 - 


AE„ 


(6) 


where cos 0 is the incident angle of the primary particle, 
and a, p and Eq are empirical parameters that need to be 
determined by a numerical simulation nH. The index/? 
describes the cutoff near the production threshold, and 
Eq is a proportionality factor applicable to the number 
of muons at the surface. In this analysis, only the pa¬ 
rameter a, describing the increase of muon multiplicity 
as a function of primary energy and mass, is relevant. 
For energies not too close to the production threshold 
Eprim/A, the relation can be simplified to: 


iV,cxA'-“.ET^ (7) 

For deep underground detectors, E^ ^in corresponds 
to the threshold energy for muons penetrating the sur¬ 
rounding material. In the case of IceCube, this corre¬ 
sponds to about 400 GeV for vertical showers, increas¬ 
ing exponentially as a function of sec flzen- 

Equationj^implies that the distribution of muon ener¬ 
gies within a shower is independent of type and energy 
of the primary nucleus, except at the very highest end 
of the spectrum, as demonstrated in Fig. The to¬ 
tal energy of the muon bundle, as well as its energy loss 
per unit track length, is therefore in good approximation 
simply proportional to the muon multiplicity. After ex¬ 
cluding rare events where the muon energy deposition is 
dominated by exceptional catastrophic losses, the muon 
multiplicity can therefore be measured simply from the 
total energy deposited in the detector. 

The experimental data can be directly related to any 
flux model expressed in terms of the parameter Emuit in¬ 
troduced in Sec. |2.3[ as long as the measured number of 
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Figure 13: Muon bundle multiplicity at closest approach 
to the center of the detector {cod) for simulated events 
with 3,000 to 4,000 registered photo-electrons. Dis¬ 
tributions are shown for trigger level and final high- 
multiplicity bundle selection. 

muons remains proportional to the overall multiplicity 
in the air shower. In the case of IceCube, the corre¬ 
sponding threshold for iron nuclei lies at about 1 PeV. 
For lower primary energies. Equation [T] is not applica¬ 
ble, and the multiplicity distribution can only be used 
for model testing, as in Sec. |4.3| 

6.2. Event Selection 

High-multiplicity bundles account for the dominant 
part of bright events in IceCube. The goal of quality 
selection is therefore not the isolation of a rare “signal”, 
but the reduction of tails and improvement in resolution. 
The criteria for the high-multiplicity bundle sample are 
shown in Table |4] 

Figure shows the true simulation-derived number 
of muons at closest approach to the center of the de¬ 
tector for events with a fixed total number of registered 
photo-electrons. On the right hand side of the distri¬ 
bution, the selection criteria eliminate very energetic 
tracks that pass through an edge or just outside the de¬ 
tector. On the left, the tail of low-multiplicity tracks 
containing high-energy muons, which are bright mainly 
because of exceptional catastrophic losses, is reduced. 

6.3. Derivation of Experimental Measurement 

The relation between the scaled parameter Emuit and 
the actual muon multiplicity in a specific detector A^^,det 
can be expressed as 

^mult = ^scale(C0S 9) ■ (8) 



Figure 14: Top: Relation between number of muons at 
closest approach to the center of the detector and to¬ 
tal energy loss of muon bundle within detector volume. 
Bottom: Total muon energy loss vs. sum of muon ener¬ 
gies at entry into detector volume. Data samples corre¬ 
spond to CORSIKA simulation after application of bun¬ 
dle selection quality criteria. The black curve represents 
a profile of the colored histogram. The error bars indi¬ 
cate the spread of the value. 


where gscaie(cos 6) is a simulation-derived function ac¬ 
counting for angular dependence of muon production 
and absorption in the surrounding material. The effects 
of local atmospheric conditions and selection efficiency 
are accounted for by a separate acceptance correction 
term. 

For the experimental measurement of the parameter 
£^muiti it is first necessary to derive expressions for the 
terms on the right hand side of of Eq. The resulting 
parameter can then be related to the analytical form of 
the bundle multiplicity spectrum by spectral unfolding. 

A numerical value of 0.79 + 0.02 for the parameter 
a was determined by fitting a power law function to 
the relation between primary energy and muon multi¬ 
plicity. The difference to the original description IflTl . 
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Selection 

Events (xlO*’) 

Rate [s 

Comment 

Effect 

All Qiot > 1,000 p.e. 

29.10 

1.075 

Base Sample (79-String Configuration) 

n/a 

cos 02en > 0.3 

28.54 

1.054 

Track Zenith Angle 

low Np 

Ldii >600 m 

24.09 

0.890 

Track Length 


high Np 

^maxlQtot ^ 0.3 

20.66 

0.763 

Brightness dominated by single DOM 

low Np 

^mpe,cod 425 m 

18.22 

0.673 

Closest approach to center of detector 

high Np 

dE/dx peak/median < 8 

12.34 

0.456 

See 

Appendix A 


low Np 


Table 4: Selection criteria and passing rates for muon multiplicity measurement. The effect of each parameter corre¬ 


sponds to a reduction at either low or high end of the distribution shown in Fig. 13 




Figure 15; Resolution of muon multiplicity estimators 
based on four different energy reconstructions. The 
analysis threshold of 1,000 photo-electrons corresponds 
to 20-30 muons. 



Figure 16: Unfolded spectra of simulated data com¬ 
pared to analytic form of spectra for three benchmark 
models uniiia. The size of the error bars corresponds 
to the expected statistical uncertainty for one year of 
IceCube data. 


Figure 17; Ratio of multiplicity spectrum unfolded sep¬ 
arately for three zenith angle regions to all-sky result. 


which gives a somewhat surprisingly accurate estimate 
of 0.757, is likely a consequence of advances in the un¬ 
derstanding of air shower physics during the last three 
decades. Recent calculations finding a lower value for 
a are only applicable in the small region of phase space 
of A ■ E/jlEprim > 0.1, where energy threshold effects 
become dominant ||5l . 

In the analysis sample, the energy loss of muons in 
the detector is in good approximation proportional to 
the number of muons Np, and to the total energy of 
muons contained in the bundle, as illustrated in Fig. 
14 An experimental observable corresponding to the 


muon multiplicity can therefore be constructed through 
a parametrization of the detector response based on 
a Monte-Carlo simulation, in the simplest case using 
the proportionality between energy deposition and total 
amount of registered photo-electrons described in Sec. 


4.1 To reduce biases and take advantage of the oppor¬ 


tunity to investigate systematic effects, the procedure 
was performed for four different muon energy estima¬ 
tors. These are: 
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• The total event charge Qtot^ measured in photo- 


















































































Figure 18; Unfolded and acceptance-corrected experimental spectrum of rescaled muon bundle multiplcity param¬ 
eter Ftnuit- The influence of systematic uncertainties listed in Table is shown separately for bin-wise fluctuations 
{uncorrelated, left) and overall scaling {correlated, right). 


Source 

Type 

Variation 

Effect 

Comment 

Composition 
Energy Estimator 
Angular Acceptance 
Light Yield 

Ice Optical 
Hadronic Model 
Seasonal Variations 
Muon Energy Loss 

uncorrelated 

uncorrelated 

uncorrelated 

correlated 

correlated 

correlated 

correlated 

correlated 

Fe, protons 

4 discrete values 

3 zenith regions 
±10% 

10% Scattering, Absorption 
discrete 

Summer vs. Winter 
Theoretical uncertainty 1681 

variable 

variable 

±10% Flux Scaling 
±13% Energy Shift 
±25% Flux Scaling 
±10% Flux Scaling 
±5% Flux Scaling 
±1% 

Residual bias near threshold 
Derived from data 

Estimated from data 

Composite Scalar Factor 

Global variations around default model 
EPOS/QGSJET/SIBYLL 
Estimated from data 

Ofhcial IceCube Value 


Table 5: Summary of systematic uncertainties affecting the result of the bundle multiplicity spectrum measurement. 



6 6.5 7 7.5 8 8.5 9 9.5 10 

log,„E„„„/GeV 



Figure 19: Interpretation of muon multiplicity spectrum by comparison to specific cosmic ray models IIT3ll^ (left), 
and by relation to all-particle flux measurements from IceCube IIZTl and other detectors ll85l[86l (right). Note that an 
exact translation to average logarithmic mass is not possible. 


electrons. Charge registered by DeepCore was ex¬ 
cluded to avoid biases due to closer DOM spacing 
and higher PMT efficiency in the sub-array. 


The mean energy deposition calculated with the 


DDDDR method described in Appendix A 


• The likelihood-based energy estimator MuEx ll64ll . 


The resolution of the multiplicity proxies in depen¬ 
dence of the true number of muons at closest approach 


• The truncated mean of the muon energy loss 
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to the center of the detector is shown in Figure [fS] Ex¬ 
cept for the raw total number of photo-electrons, all es¬ 
timators perform in a remarkably similar way in sim¬ 
ulation. The presence of individual outliers illustrates 
the motivation to use more than one method to ensure 
stability of the result. The angular-dependent scaling 
function gscaie(cos 0zen) was parametrized based on sim¬ 
ulated data. 

Using the RooUnfold algorithm Il84ll . a spectral un¬ 
folding was applied to the measured distribution of 
E,„uit. The differential flux was then related to the un¬ 
folded and histogrammed experimental data as; 




— c(A£’bin5 ^exp) ' ^(Emult) 


AiVe, 


^ logjQ Emult 


(9) 


Here the proportionality constant c accounts for the 
effective livetime of the data sample and the bin size of 
the histogram. The detector acceptance //(Emuit), whose 
exact form depends on atmospheric conditions, needs to 
be derived from simulation. 

The approach can be verified by a full-circle test, as 


shown in Fig. 16 Each of the benchmark models, cho¬ 
sen to reflect extreme assumptions about the behavior of 
the cosmic ray flux, can be reproduced by applying the 
analysis procedure to simulated data. 

6.4. Result 

Systematic uncertainties applying to the experimen¬ 
tal measurement are summarized in Table |5] The cate¬ 
gorization by type corresponds to bin-wise fluctuations 
{uncorrelated) and overall scaling effects {correlated). 

Of special interest is the angular variation, which 
dominates the total bin-wise uncertainty over a wide 
range. The effect is illustrated in Fig. [m Splitting 
the data according to the reconstructed zenith angle 
into three separate event samples results in spectra that 
are similar in shape, but whose absolute normalization 
varies within a band of approximately +10%. As the 
difference appears to be not uniform, it has been conser¬ 
vatively assumed to lead to uncorrelated bin-wise varia¬ 
tions in the all-sky spectrum. Notwithstanding, magni¬ 
tude and direction are similar to the unexplained effect 
described in Section suggesting a possible common 
underlying cause. The final result, after successive ad¬ 
dition of systematic error bands in quadrature, is shown 


in Fig. 18 


Since the muon multiplicity is not a fundamental pa¬ 
rameter of the cosmic ray flux, it is important to find 
an appropriate way for its interpretation. Two possibil¬ 
ities are illustrated in Fig. The first is by express¬ 
ing cosmic ray flux models in terms of Emuit through 


application of Eq. [T] Experimental result and predic¬ 
tion can then be directly related. The second is to trans¬ 
late the multiplicity distribution to an energy spectrum 
under a particular hypothesis for the elemental compo¬ 
sition. By default, the scaling of Fmuit corresponds to 
iron, but changing it to any other primary nucleus type 
is straightforward. 

The result can then be overlayed by independent cos¬ 
mic ray flux measurements. An unambiguous derivation 
of the average mass as a function of the primary energy 
is not possible due to the degeneracy between mass and 
energy in the multiplicity measurement. However, the 
qualitative variation of composition with energy is con¬ 
sistent with a gradual change towards heavier elements 
in the range between the knee and 100 PeV. If the cur¬ 
rent description of muon production in air showers is 
correct, and the external measurements are reliable, a 
purely protonic flux would be strongly disfavored up un¬ 
til the ankle region. 


7. High-Energy Muons 


7.1. Principle 

The presence of a single exceptionally strong catas¬ 
trophic loss can be used both for tagging high-energy 
muons and to estimate their energy. The first part is 
obvious; An individual particle shower along a track 
can only have been caused by a parent of the same en¬ 
ergy or above. Simulated data indicate that instances in 
which two or more muons in the same bundle suffer a 
catastrophic loss simultaneously in a way that is indis- 
tiguishable in the energy reconstruction are exceedingly 
rare. 

The quantification is based on the close relation be¬ 
tween the energy of the catastrophic loss used to identify 
the event and that of the leading muon, a consequence 
of the steeply falling spectrum. Once the muon energy 
at the point of entry into the detector volume has been 
determined, the most likely energy at the surface of the 
ice can be estimated by taking into account the zenith 


angle, as illustrated in Fig. 20 This method was de¬ 


veloped specifically for the purpose of measuring the 
energy spectrum of atmospheric muons. As shown in 
Fig. [T^ the leading particle typically only accounts for 
a limited fraction of the total event energy, and the ap¬ 
plication of energy measurement techniques optimized 
for single neutrino-induced muon tracks could lead to 
substantial biases in the case of a large accompanying 
bundle. 

Higher-order corrections are necessary to account for 
correlations and the effect of variations in the distance to 
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Figure 20; Top; Relation between most energetic sin¬ 
gle energy loss and leading muon energy within the Ice- 
Cube detector volume. Middle; Distribution of true en¬ 
ergy parameters for two slices in top histogram. Bot¬ 
tom; Fraction of muon surface energy remaining at 
point of entry into detector volume in dependence of 
zenith angle. Figures are based on simulated events with 
primary flux weighted to ^ power law spectrum and 
correspond to hnal analysis sample before application 
of minimum shower energy criterion. The black curves 
represent mean and spread of the distribution. 


the surface due to the vertical extension of the detector. 
All relations in this study were based on parametriza- 
tions using simulated events. A full multi-dimensional 
unfolding would be preferable, but requires a substantial 
increase in simulation statistics. 

7.2. Event Selection 




Figure 21; Example for peak to median energy loss ra¬ 
tio in high-energy muon candidate event found in ex¬ 
perimental data. Top; Reconstructed differential en¬ 
ergy loss in dependence of distance to surface, measured 
along the reconstructed track. Details of the method 
are described in Appendix A Bottom; Image of the 
event. The volume of each sphere is proportional to the 
signal registered by a given DOM. The color scheme 
corresponds to the arrival time of the hrst photon (red; 
earliest, blue; latest). Reconstructed event parame¬ 
ters are; TeV, = 1.03!° “ PeV, 

0zen = 45.1 + 0.2° 


The selection of muon events with exceptional 
stochastic energy losses is primarily based on recon¬ 
structing the differential energy deposition and selecting 
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tracks according to the ratio of peak to median energy 
loss as illustrated in Fig. 21 All other criteria are an¬ 
cillary, and are only applied to minimize a possible con¬ 
tribution from misreconstructed tracks. An overview of 
the selection is given in Table 


A special case is the exclusion of events with a recon¬ 
structed shower energy of less than 5 TeV. This require¬ 
ment was added to reduce uncertainties in the thresh¬ 
old region, which may not be well described by current 
understanding of systematic detector effects. The rea¬ 
son to choose a value of 5 TeV is that a typical electro¬ 
magnetic shower of that energy will produce a signal of 
about 1,000 photo-electrons, coinciding with the base 
sample selection. 


logio£casc,reco/GeV = 1.6888 ■ eO^l^'og.o Eea..„„/GeV 

( 10 ) 

In the energy region between 5 TeV and 1 PeV, the 
difference between raw and final value is smaller than 
0.1 in logjo E. 

The stochastic energy loss ficasc.reco was then used to 
estimate the most likely energy of the leading muon at 
the surface in dependence of zenith angle 0zen and 
slant depth t/siant = Zvert/ cos 6>zen, where Zvert is the ver¬ 
tical distance to the surface at the point of closest ap¬ 
proach to the center of the detector. 


7.3. Energy Estimator Construction 



Figure 22: Relation between reconstructed and true sur¬ 
face energy for simulated atmospheric muon data before 
excluding events with reconstructed shower energy of 
less than 5 TeV. The primary particle flux in the simula¬ 
tion was weighted according to a power law of the form 
. Also shown are mean and spread of the distribu¬ 
tion. 


The energy reconstruction is based on the determin¬ 


istic reconstruction method discussed in Appendix A 


which was designed specifically for this purpose. Sub¬ 
sequently developed likelihood methods Il64l were eval¬ 
uated, but gave no improvement in resolution while in¬ 
troducing a tail of substantially overestimated energies. 

In the first step, the energy ficasc.reco of the strongest 
loss (“cascade”) along the track was determined. The 
exact value is almost identical to the raw reconstructed 
energy ficasc.raw from the DDDDR algorithm, except for 
a minor correction factor of the form: 



Figure 23: All-Sky surface flux predictions BTI for 
three different cosmic ray models and spectrum ex¬ 
tracted from full IceCube detector simulation with same 
primary weight. The error bars on the measured spec¬ 
trum are the consequence of limited statistics. 


The parametrized form of the measured muon surface 
energy is: 


logio£;“ eo/GeV = 0.554 -h 0.884- 

' ^casc,reco/GeV) + ^zen? ^slant)) 

( 11 ) 


where /corr(cos 0zen,'^siant) Is a fifth-oi'der polynomial. 
This relation represents a purely empirical parametriza- 
tion based on the interpolation of detector-specific sim¬ 
ulated data. 


The relation between the experimental muon surface 
energy estimator defined in Eq. 11 and the true energy 
of the leading muon at the surface is shown in Fig. 

It is important to note that the definition is only valid for 
spectra reasonably close to that used in the construction. 
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Quality Level 

Events (xlO^) 

Rate [i- *] 

Comment 

All 2tot > 1 > OOOp.e. 

38.28 

1.334 

Base Sample (86-String Configuration) 

cos >0.1 

37.99 

1.324 

Track zenith angle 

^msx/Qtot 0.5 

34.46 

1.201 

Brightness dominated by single DOM 

Ljir > 800m 

27.55 

0.960 

Track length in detector 

^DOM.lSOm > 40 

24.71 

0.861 

Stochastic loss containment 

peak/median dEjdx > 10 

2.795 

0.0974 

Exceptional energy loss along track 

median dEjdx > 0.2GeV/m 

2.769 

0.0965 

Exclude dim tracks 

Ecasc > 5 TeV 

0.769 

0.0268 

Exclude threshold region 


Table 6; High-energy muon selection criteria and passing rates. 




Figure 24: Experimentally measured spectrum of high-energy muons using one year of IceCube data. Left: All-Sky 
flux with bin-by-bin and correlated error margin. Right: All-Sky spectrum compared to flux above and below 60 
degrees (cos 0zen = 0.5). Only bin-wise errors are shown. Between 15 TeV and 1.5 PeV, the all-sky spectrum is 
consistent with a power law of index = -3.78, illustrated by the dashed line. 


Source 

Type 

Variation 

Effect 

Comment 

Composition 
Angular Acceptance 
DOM Efficiency 
Optical Ice 
Seasonal Variations 
Muon Energy Loss 

uncorrelated 

uncorrelated 

correlated 

correlated 

correlated 

correlated 

Fe, protons 

0.2 • (cos 0zen “ 0.5) 
±10% 

10% Scattering, Absorption 
Summer vs. Winter 
Theoretical uncertainty l68l 

variable 

See Text 

±10% Energy Shift 
±10% Energy Shift 
±5% Flux Scaling 
±1% 

Negligible Above 25 TeV 
Unknown Cause 
Effective light yield 
Clobal variations 
Prompt Invariant 
Official IceCube Value 


Table 7: Systematic uncertainties in the calculation of the high-energy muon energy spectrum. 


7.4. Energy Spectrum 


The final muon energy spectrum was calculated by 
dividing the histogrammed number of measured events 
T^data by a generic prediction from a full detector sim¬ 
ulation A^detMC, and then multiplying the ratio with the 
corresponding flux OsmfMC at the surface. Specifically, 
IceCube detector simulation and external surface data 
set ED were weighted according to a power law of the 
formE-^-’: 


(i® 


/j,exp 


dEu 


data 
A psurf 

^^/j,reco 


AA^detMC2./ 

A psurf 

^^/j,reco / 




'f/,surfMC2.7 


dE 


surf 


( 12 ) 


Figure 23 demonstrates the validity of the analysis 


procedure, and the robustness of the energy estima¬ 
tor constrTiction against small spectral variations. The 
surface flux for different primary model assumptions 
can be extracted accurately from simulated experimen¬ 
tal data. While a full unfolding would be preferable, the 
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currently available simulated data statistics do not allow 
for the implementation of such a procedure. 

In the derivation of the experimental result, the sys¬ 
tematic uncertainties listed in Table[7]were applied. The 
classification according to correlation is the same as 
in Section 6.4 Except for a small effect due to pri¬ 
mary composition near threshold, all experimental un¬ 
certainties lead to correlated errors. A special case 
is the angular acceptance. In light of the low-energy 
muon and multiplicity spectrum studies described in 
Sec. |4.3[ it is necessary to take into account the possi¬ 
bility of an unidentified error source distorting the dis¬ 
tribution. This was done by calculating the energy spec¬ 
trum once for the default angular acceptance and once 
with simulated events re-weighted by an additional fac¬ 
tor Wcorr = cr- (cos 0zen -0.5), where a corresponds to an 
ad-hoc linear correction parameter. The value a - 0.2, 
corresponding to the variation of + 10% seen in the other 
analyses, reflects the assumption that the effect is inde¬ 
pendent of the event sample. 

The experimentally measured muon energy spectrum 
is shown in Fig. Distortion due to possible angu¬ 
lar effects are small compared to the statistical uncer¬ 
tainty. Within the present accuracy, the average all-sky 
flux above 15 TeV can be approximated by a simple 
power law: 


—a = 1.06+SS X lO^'^s^^cm^^srad^^TeV^^ 

/ p X-3.78±0.02(ifaf.)±0.03(j>'sf.) 

(lolv) 

The translation to a vertical flux as commonly used in 
the literature is not trivial, since the angular dependence 
of the contribution from prompt hadron decays is differ¬ 
ent from that of light mesons, and its magnitude a priori 
unknown. 

The almost featureless shape of the measured spec¬ 
trum might appear as a striking contradiction to the 
naive expectation of seeing a clear signature of the sharp 
cutoff of the primary nucleon spectrum at the knee. 
However, closer examination reveals that this is very 
likely a simple coincidence resulting from the fact the 
the prompt contribution approximately compensates for 
the effect of the knee if the flux is averaged over the 
whole sky. 

Calculating the spectra separately for angles above 
and below 60 degrees from zenith shows the expected 
increase of the muon flux toward the horizon. Beyond 
approximately 300 TeV, the two curves appear to con¬ 
verge, consistent with the emergence of an isotropic 


prompt component. A quantitative discussion of the an¬ 
gular distribution is given in the following section. 

The final all-sky spectrum was then fitted to a combi¬ 
nation of “conventional” light meson and prompt com¬ 
ponents, with a Gaussian prior of Ay = 0.1 applied to 
the spectral index. The result in the case of H3a and 
GST-GF models is illustrated in Fig.|^ The difference 
between the two measurements is due to the presence 
of a spectral component in the GST-GF model with a 
power-law index of -2.3 to -2.4 compared to about -2.6 
in H3a. Even though the exponential cutoff energy of 4 
PeV is identical in both cases, the influence of the steep¬ 
ening at the knee is effectively reduced in the harder 
spectrum. 


The best fit values for the prompt contribution are 
listed in the second column of Table |8] relative to the 
ERS flux 1531. Note that unlike the theoretical predic¬ 
tion, which applies specifically to neutrinos from charm, 
the experimental result presented here is the sum of 
heavy quark and light vector meson decays. A detailed 


discussion can be found in Appendix B 


Since only the energy spectrum is used here, the par¬ 
tial degeneracy between the behavior of the all-nucleon 
flux at the knee and the prompt contribution is pre¬ 
served. Consequently, the magnitude of the prompt 
component strongly depends on the primary model. Ex¬ 
cept for the proposal by Zatsepin and Sokolskaya llTSl . 
each of the flux assumptions can be reconciled with the 
data without a major spectral adjustement. 


7.5. Angular Distribution 


The ambiguity between nucleon flux and prompt con¬ 
tribution can be resolved by the addition of angular in¬ 
formation. Figure]^ shows the best fit results from the 
previous section compared to data separately for angles 
above and below 60 degrees from zenith. While neither 
of the two models shown here is obviously favored, it 
is clear that a substantial prompt contribution is needed 
in either case to explain the difference between the two 
regions. 


A quantitative treatment can be derived from the dif¬ 
ferent behavior of light meson and prompt components. 
The prompt flux is isotropic, whereas the contribution 
from light meson decays is in good approximation in¬ 
versely proportional to cos 0zen 61. Using the prompt 
flux description derived in Appendix B the experimen¬ 
tally measured fraction of prompt muons as a function 
of muon energy and zenith angle is: 
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Figure 25; All-sky muon energy spectrum and predictions based on H3a (left) and Global Fit model (right) IT3l . Best 
fit parameters are listed Table 


CR Model 

Best Pit (ERS) 

A^/dof 

Icr Interval (90% CL) 

Pull (Ay) 

(^■(^Prompt > 0) 

GST-Global Eit ^311 

2.14 

7.96/9 

1.27-3.35 (0.77-4.30) 

0.01 

2.64 

H3a lia 

4.75 

9.09/9 

3.17 -7.16 (2.33 -9.34) 

-0.03 

3.97 

Zats.-Sok. 135] 

6.23 

13.98/9 

4.55 - 8.70 (3.59 - 10.68) 

-0.23 

5.24 

PG Constant Ay ||33l 

0.94 

9.07/9 

0.36- 1.63 (<2.15) 

0.03 

1.52 

PG Rigidity 1331 

6.97 

5.86/9 

4.73 - 10.61 (3.53 - 13.83) 

-0.06 

4.35 


Table 8: Result of model-dependent fit to all-sky muon energy spectrum. Note that for muons, the prompt flux 
is expected to include a substantial contribution from electromagnetic decays of light vector mesons, which is not 
present in neutrino spectra ll56l . 


/prompt(^/j? cos G) — 


^prompt (^/r AOS 6 ^ 

<l>total(£';j, COS 0) 


- 1 + 


£ 1/2 ■ COS 0 
■ /corr(£/r) 


(14) 


In this approximation, the prompt contribution is de¬ 
scribed independent of the muon flux ^^(ii^). The 
repartition between the two components at a given en¬ 
ergy can therefore be measured from the angular distri¬ 
bution alone. The effect of higher order terms, such as 
departure of the angular distribution from a pure sec 0zen 
dependence due to the curvature of the Earth and devia¬ 
tions of the nucleon spectrum from a simple power law, 
have been estimated as less than 10% using a full DP- 
MJET ||87l simulation of the prompt component. 

In this study, the measurement of the prompt flux was 
based on splitting the event sample into two separate 
sets according to the reconstructed zenith angle. The 
ratios between experimental data and Monte-Carlo sim¬ 
ulation were then combined into a single parameter de¬ 
fined as: 



Eigure 27: Ratio parameter rhor.ven expressing deviation 
of angular distribution from purely conventional flux for 
various prompt levels in simulation. The size of the er¬ 
ror bars corresponds to the statistical uncertainty due to 
limited availability of simulated data. 
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Figure 26; Horizontal and vertical muon energy spectra compared to prediction using best fit values to all-sky spectrum 
as listed in Table[^ Top row: vertical (0-60 degrees from zenith), bottom row; horizontal (60-84 degrees from zenith). 
Left; H3a, Right; Global Fit Model. 


^hor.vert “ 


A^,z,data(6'zen > 60°) /A/;,,data(6'zen < 60°) 


-^/i,Mc(®zen > 60°) \A^^,Mc(0zen < 60°) 


-1 


(15) 

The variation as a function of muon energy is illus¬ 
trated in Fig. 27 where N^^mc represents the purely 
conventional flux, and A)j,data is derived from simula¬ 
tion weighted according to two assumptions about the 
prompt flux level. Using two discrete samples is not the 
most statistically powerful way to exploit the angular 
information, but minimizes fluctuations resulting from 
limited simulation availability. 

The experimental result is shown in Fig. The 
best estimate for the prompt flux is significantly higher 
than the theoretical prediction, but well within the mar¬ 
gin permitted by the model-dependent fits to the energy 
spectrum discussed in the previous section. 

Given the presence of an unknown systematic error in 
the low-level and high-multiplicity atmospheric muon 
samples as described in Sec. |4.3| it is necessary to take 
into account the possibility that the angular distribution 
might be distorted. As the source of the effect is still 
unknown, the only choice is to evaluate the influence on 
the measurement by applying a generic correction term. 

Figurej^shows the consequence of re-weighting the 



Figure 28: Best angular prompt fit using default as¬ 
sumptions about systematic uncertainties. Expressed in 
multiples of the ERS flux 1531 . the result is 4.9 + 0.9, 
with ;if2/dof=20.0/15. 
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Figure 29: Top: Two-dimensional probability distribu¬ 
tion function of angular prompt fit results in the pres¬ 
ence of an ad-hoc correction term as described in Sec¬ 
tion 7.4 The y-axis corresponds to the angular adjust¬ 


ment parameter a. Bottom: Result for best overall fit 
with;t^Vdof=14.9/15, located at (2.41; 0.18). 


simulated data by a linear term of the form 1 H- a ■ 
(cos 0zen)- The two-dimensional distribution demon¬ 
strates that an imbalance between horizontal and ver¬ 
tical tracks with a magnitude of 18% describes the data 
best. This value is suggestively close to the distortions 
observed in Sec. |^and[^ although the limited statistical 
significance does not permit a firm conclusion. 

7.6. Discussion 

A definite measurement of the prompt flux is not yet 
possible. Depending on which assumption is chosen for 
the systematic error, the final result varies considerably. 
Figure shows the significance levels for default as¬ 
sumption and full marginalization over the linear cor¬ 
rection factor. Best fit values and confidence intervals 
for each case are listed in Tabled 

At present, the best neutrino-derived limit for the at¬ 
mospheric prompt flux is 2.11 ERS at 90% confidence 


Figure 30: Significance of prompt flux measurement 
based on angular information. The individual curves 
correspond to different assumptions about systematic 
effects as described in the text. Also shown is the hy¬ 
pothetical result which could be achieved with one year 
of experimental data given unlimited availability of sim¬ 
ulated events, assuming a best fit value of 1.8 ERS con¬ 
sistent with theoretical predictions for inclusive prompt 
muon flux. 


level 1881 . This result was derived by a likelihood fit 
combining four independent measurements from Ice- 
Cube, and includes both track-like (v^ charged current) 
and shower-like (Vg and charged current, all-flavor 
neutral current) neutrino event topologies. Eor com¬ 
parisons it is important to keep in mind that the at¬ 
mospheric muon measurement result represents the in¬ 
clusive prompt flux, potentially including a substantial 
contribution from electromagnetic decays of unflavored 
vector mesons l5^ . It is also worth noting that recent 
studies show that the uncertainty of theoretical models 
for atmospheric lepton production in charm decays are 
larger than previously assumed l89l . 

None of the model fluxes selected for the fit to the 
muon energy spectrum requires a prompt flux in dis¬ 
agreement with the neutrino measurement, with the ex¬ 
ception of the proposal by Zatsepin and Sokolskaya. 
The rigidity-dependent poly-gonato model lacks an ex- 
tragalactic component whose inclusion would lead to a 
higher nucleon flux and therefore a lower estimate for 
the prompt contribution. 

The result based on the angular distribution alone is 
almost independent of the nucleon flux and would even 
at the present stage be statistically powerful enough 
to constrain competing primary nucleon flux models 
around the knee. Unfortunately this possibility is pre¬ 
cluded by the likely presence of an unidentified sys¬ 
tematic error source. Both uncorrected and ad-hoc cor- 
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Sample 

Best Fit (ERS) 

Icr Interval (90% CL) 

0'(<t)prompt > 0) 

Uncorrected 

4.93 

4.05-5.87 (3.55-6.56) 

9.43 

Marginalized Ang. Corr. 

3.19 

1.64-5.48 (0.98-7.26) 

3.46 


Table 9; Result of Angular Prompt Fit. 


rected measurements could be reconciled with different 
predictions based on data from air shower arrays, no¬ 
tably the H3a and Global Fit models na. At present, 
the angular measurement is also fully consistent with 
constraints derived from neutrino data. 

8. Conclusion and Outlook 

The influence of cosmic rays on IceCube data is sig¬ 
nificant and varied. Given the presence of several en¬ 
ergy regions where external measurements by direct de¬ 
tection or air shower arrays are sparse, it is necessary 
to develop a comprehensive picture including neutrinos, 
muons and surface measurements. Atmospheric muons 
play a privileged role, as they cover the largest energy 
range and provide the highest statistics. A consistent 
description of all experimental results will be an impor¬ 
tant contribution for the understanding of cosmic rays 
in general. 

The studies presented in this paper have outlined the 
opportunities to extract meaningful results from atmo¬ 
spheric muon data in a large-volume underground par¬ 
ticle detector. Once systematic effects are fully under¬ 
stood and controlled, it will be possible to measure the 
muon energy spectrum from 1 TeV to beyond 1 PeV 
by combining measurements based on angular distribu¬ 
tion and catastrophic losses. Agreement between the 
two methods can then be verified in the overlap region 
around 10-20 TeV. 

There is a strong indication for the presence of a com¬ 
ponent from prompt hadron decays in the muon energy 
specrum, with best fit values generally falling at the 
higher side of theoretical predictions. In the future, it 
will be possible for the IceCube detector to precisely 
measure the prompt contribution and to constrain the 
all-nucleon primary flux before and around the knee. 
With more data accumulating, independent verification 
of the prompt measurement based on seasonal variations 
of the muon flux 1901 will soon become feasible as well. 

The muon multiplicity spectrum provides access to 
the cosmic ray energy region beyond the knee. Even 
though a direct translation of the result to primary en¬ 
ergy and average mass is impossible, combination with 
results from surface detectors or comparisons to model 


predictions provide valuable insights. In coming years, 
the measurement can be extended further into the tran¬ 
sition region around the ankle. A possible contribution 
from heavy elements to the cosmic ray flux at EeV en¬ 
ergies should then be discernible. 

An important goal of this study was to verify the cur¬ 
rent understanding of systematic uncertainties. An un¬ 
explained effect was demonstrated using low-level data, 
and appears to be present in the other analysis samples 
as well. In order to improve the quality of future at¬ 
mospheric muon measurements with IceCube, it will be 
essential to determine whether the observed discrepancy 
requires better understanding of the detector, or of the 
production mechanisms of muons in air showers. 

Comparisons with measurements from the upcoming 
water-based KM3NeT detector El] will be invaluable 
to decide whether the inconsistencies seen in IceCube 
data are due to the particular detector setup, or represent 
unexplained physics effects. 
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Appendix A. Data-Derived Deterministic Differ¬ 
ential Deposition Reconstruction 
(DDDDR) 

Appendix A.l. Concept 

The energy deposition of muons at TeV energies 
passing through matter is not continuous and uniform, 
but primarily a series of discrete catastrophic losses. 
In order to exploit the information contained in the 
stochasticity of muon events, it is necessary to recon¬ 
struct the differential energy loss along their tracks. The 
study presented in this paper requires a robust method 
for identification and energy measurement of major 
stochastic losses. Its principle is to use muon bundles 
in experimental data to characterize photon propagation 
in the detector and apply the result to the construction 
of a deterministic energy estimator. 



Figure A.l: Sketch of light attenuation around muon 
track in ice. 


Figure A.l shows a sketch of the photon intensity 
distribution around the reconstructed track of a muon 
bundle. In the ideal case of a perfectly transparent 
homogeneous medium and a precisely defined infinite 
one-dimensional track of arbitrarily high brightness, the 
light intensity would fall off as 1 /t/jp, where the impact 
parameter dip is defined as the perpendicular distance 
to the track. Assuming the measured charge ^dom in 
a given DOM to be proportional to the light density, 
and the emitted number of photons Aphot to be propor¬ 
tional to the energy deposition the relation be¬ 
tween muon energy deposition and measurement then 
takes the form: 


AEf^lAx ~ Aphot ~ ?DOM ■ (A.l) 

In reality, scattering and absorption in the detector 
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medium require the addition of an exponential attenu¬ 
ation term exp(-t/ip//latt): 


set fully covered by the online event filters. For each 
DOM within a given vertical depth range, the quantity 


A^phot ^DOM ■ dw ■ expidi^l (A.2) 

where the attenuation length datt depends on the lo¬ 
cal optical properties in a given part of the detector. 
Approximating the structure of individual ice layers as 
purely horizontal, datt is simply a function of the vertical 
depth Zvert- 
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Figure A.2: Top: Lateral attenuation of photon intensity 
along muon bundle tracks in experimental data. The 
vertical depth ranges, corresponding to DOM position 
relative to the center of the detector 1949 m below the 
surface, were chosen to illustrate the strongly varying 
optical properties of the ice. Bottom: Effective attenua¬ 
tion parameter Tatt derived from exponential fit to the 
data distribution. Experimental values are compared 
to Monte-Carlo simulation using reconstructed and true 
track parameters for calculation of the impact parameter 

diy. 


The validity of this hypothesis is demonstrated in 
Eig. |A.2| A sample of bright downgoing tracks with 
2tot > 1, OOOpe was selected to obtain an unbiased data 


^phot,ideal ~ ^DOM ' ^DOM ' ^IP (-^-3) 

is calculated, corresponding to the photon yield ad¬ 
justed for the distance from the track and relative quan¬ 
tum efficiency edom of the PMT, which is 1 in standard 
and about 1.35 for high-efficiency DeepCore DOMs. 
The curves are averaged over the entire event sample 
and include DOMs that did not register a signal. The 
solid lines shows the result of a fit to the function 


/(c/ip) = c ■ e.rp(-c/ip/datt) (A.4) 

with the effective attenuation length Tatt and the data 
sample-dependent normalization constant c as free fit 
parameters. Exponential attenuation as a function of 
the impact parameter is a valid assumption over a wide 
range, breaking down only for very close distances and 
in the layer with high dust concentration at Zven ~ 
-100 m, where the vertical gradient of the optical ice 
properties is exceptionally steep. 

The experimental result is well reproduced by the 
simulation, as illustrated in the lower plot. The very 
small difference between the curves using true and re¬ 
constructed track parameters means that track recon¬ 
struction inaccuracies can be neglected. 

Appendix A.2. Construction of Energy Observable 

Once the effective attenuation length has been deter¬ 
mined, it can be used to construct a simple differential 
energy loss parameter. Eor each DOM within a given 
distance from the reconstructed track, an approximation 
for the photon yield corrected for PMT efhciency and 
ice attenuation can be calculated. The actual differen¬ 
tial energy loss at the position of the DOM projected is 
related to the experimental observable by: 


/ 

\ dx /, 


%OM ■ 9'DOM' 


u, 

1 ^IP " ^ 


d\^ < do 
d\p > do 


(A.5) 


where /scale - 0.020GeV ■ (p.e is a simple scal¬ 
ing factor that can be derived from a Monte Carlo simu¬ 
lation and dQiz) - 19mH-0.01 -z expresses the mild depth 
dependence of the point of transition from flat to expo¬ 
nential behavior. The vertical coordinate z is measured 
from the center of the detector at 1949 meters below the 
surface. 
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To account for fluctuations affecting individual mea¬ 
surements and DOMs that did not register a signal, the 
track is subdivided into longitudinal bins with a width of 
50 meters, over which the measured parameter is aver¬ 
aged. The lateral limit for the inclusion of DOMs can be 
adjusted to find a compromise between sufficient statis¬ 
tics and adequate longitudinal resolution. The principle 
is illustrated in Fig. |A.3| 

Note that the exact value of dE/dx is only calcu¬ 
lated for demonstration purposes and should be con¬ 
sidered approximate. The measured observables, like 
any energy-dependent observable, are in practical ap¬ 
plications directly related to physical parameters such 
as shower energy and muon multiplicity, where the ex¬ 
act conversion depends on the spectrum of the data dis¬ 
tribution. 

The energy of the strongest stochastic loss in the 
event could be derived immediately from the highest bin 
value in the profile. However, this estimate is often im¬ 
precise. Better results can be achieved by a dedicated 
reconstruction for the individual loss energy. The origin 
of the shower is assumed to coincide with the position 
of the DOM with the highest dEjdx value projected on 
the track. Its energy is then calculated in a similar way 
as for the track, except that the photon emission is as¬ 
sumed to be point-like and isotropic. Instead of falling 
off linearly, the light intensity falls off quadratically as a 
function of distance, and the energy estimate becomes: 

Floss,reco ~ ^DOM ' ^DOM' 

, jrl, Hoss < ro (A.6) 

noss>ro 

The shower energy can then be determined by calcu¬ 
lating the mean of the values for the individual DOMs. 
The energy resolution for events selected by the method 
described in Section|7]is shown in Fig. |A.4| 

Appendix B. Prompt Flux Calculation 

Appendix B.l. Prompt Muon Flux Approximation 

The characteristics of the atmospheric muon energy 
spectrum at energies beyond 100 TeV are influenced by 
prompt hadron decays. In neutrino analyses, these can 
be taken into account by applying a simple weighting 
function to simulated data. Muons, on the other hand, 
are always part of a bundle, and in principle it would 
be necessary to generate a full air shower simulation in¬ 
cluding prompt lepton production. 

The hadronic interaction generators integrated into 
the CORSIKA simulation package as of version 7.4 are 


not adequate for a prompt muon simulation mass pro¬ 
duction. QGSJET and DPMJET lIFTI are slow, and 
charm production in QGSJET is very small compared 
to theoretical predictions. The core CORSIKA prop¬ 
agator does not handle re-interaction effects for heavy 
hadrons, which become important at energies approach¬ 
ing 10 PeV. 

A version of SIBYLL that includes charm is at the 
development stage 1921 . The updated code also takes 
into account production and decay of unflavored light 
mesons, which form an important part of the prompt 
muon flux l54l . First published simulated prompt atmo¬ 
spheric muon spectra indicate consistency with the ERS 
model for charmed mesons, and an unflavored compo¬ 
nent of approximately equal magnitude ll5^ . 

In this paper, the prompt flux is expressed in depen¬ 
dence of the “conventional’ flux from light meson de¬ 
cays. In this way it can be modeled using simulated 
events from the standard IceCube CORSIKA mass pro¬ 
duction, including detector simulation and information 
about the primary cosmic ray composition. 

Construction of the simulated prompt flux is based on 
the following assumptions: 

• The spectral index of the prompt component 
Tprompt is related to the conventional index yconv as 
Tprompt = Tconv + 1- Highei'-ordef effects, such as 
the varying cross section of charm production and 
re-interaction in the atmosphere, can be accounted 
for by a corrective term fconiEfj). 

• The prompt flux is isotropic, the conventional flux 
increases proportional to sec flzen in the analysis re¬ 
gion above cos 0zen = 0.1. Variations due to the 
curvature of the Earth 15^ are neglected. 

• The influence of changes in the nucleon spectrum 
on the prompt flux is the same as on the con¬ 
ventional flux. Based on estimates using prompt 
muons simulated with DPMJET, this assumption 
is valid within 10% for spectra with an exponential 
cutoff at the knee. 

• The contribution from light vector meson di-muon 
decays is small compared to that from heavy 
hadrons and/or has the same energy spectrum. For 
prompt muon fluxes simulated with the newest de¬ 
velopment version of SIBYLL, charm and unfla¬ 
vored spectra are almost identical in shape between 
10 TeV and 1 PeV fSh). 

The approximated prompt flux is then: 
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^/J,prompt(^/i? ^zen) ~ ^/j,conv(^/j5 ^zen) 


E/j ■ COS 0z, 

£ 1/2 


' fcon(Eu) 


(B.l) 


The relative flux normalization is expressed in terms 
of £ 1 / 2 , the crossover energy for prompt and conven¬ 
tional fluxes in vertical air showers. This parameter pro¬ 
vides a simple and intuitively clear way to express the 
magnitude of the prompt flux, and can easily be esti¬ 
mated. 

To calculate the crossover energy £ 1/2 for a specific 
prediction, it is sufficient to compare conventional muon 
simulations with a prompt flux parametrization, as il¬ 
lustrated in Fig. B.l The crossover energy can then 


be determined in a straightforward way by a fit to their 
ratio. Note that here the primary nucleon spectrum cor¬ 
responds to the naive TIG model Il93ll used in the theo¬ 
retical calculation. 

Since the full air shower simulation only needs to pro¬ 
vide an estimate for the conventional flux, this proce¬ 
dure can be repeated for any interaction model. In this 
study, as in most IceCube analyses, the prompt predic¬ 
tion is based on the calculation by Enberg, Reno and 
Sarcevic 15^ . The corresponding values are listed in 
Table iBTl 


Hadronic Model 

ERS (max) 

ERS (default) 

ERS (min) 

SIBYLL 

QGSJET-II 

QGSJET-OIc 

5.71 ±0.02 
5.62 ± 0.02 
5.65 ± 0.02 

5.82 ±0.03 
5.72 ± 0.03 
5.75 ± 0.03 

5.99 ± 0.03 
5.90 ± 0.03 
5.93 ± 0.03 


Table B.l: Vertical crossover energy logjo£i/ 2 /GeV for 
ERS flux and CORSIKA non-prompt muon simulation. 


Detailed features of a theoretical model are taken 
into account by a higher-order correction. In partic¬ 
ular, those are the increase of the prompt production 
cross section as a function of primary energy and the ap- 
pearence of re-interaction effects at energies of several 
PeV. Since the latter is negligible in the range covered 
by the study in this paper, its angular dependence was 
omitted. 

The parametrized form of the correction factor is: 


/co/'r(£/j) “ •') ' ./corr(^'^f) — 

[(3.74 - 0.461 ■ logio£^/GeV) ■ (1 -H e2 i3 i°gioEV4.9PeV^j-' 

(B.2) 


After application of the correction, simulation-based 
flux prediction and theoretical model agree well, as il¬ 


lustrated in Eig. B.2 


Appendix B.l. Translation to Neutrino Flux 

Prompt muon and neutrino fluxes are not strictly 
identical. In particular, muons can originate in elec¬ 
tromagnetic di-muon decays of vector mesons. The 
muon-derived measurement is a combination of unfla¬ 
vored and heavy quark-induced fluxes: 

®prompt,/j = ®/i,heavy + ®unflav (B.3) 

Whereas previous estimates based on theoretical cal¬ 
culations indicated an unflavored contribution of 0.3-0.4 
times the ERS flux Ea, recent numerical simulations 
result in a higher value, almost approaching the flux 
from heavy hadron decays ll56ll . 

The contribution from vector meson decays is par¬ 
tially compensated by a relative suppression of the 
muon flux with respect to neutrinos of 15-20% originat¬ 
ing in the physics of c —> i decay Il42l . here represented 
by the conversion factor The resulting neutrino flux 
is therefore: 


^prompt,V “ ^/J,v ' (^prompt,/j ^unflav) (B.4) 

An exact translation requires precise determination of 
spectrum and magnitude of the unflavored contribution 
and evaluation of the weak matrix element responsible 
for At the moment, the calculation of a reliable es¬ 
timate for the prompt atmosperic neutrino flux is pre¬ 
cluded by the substantial uncertainties on the experi¬ 
mental measurement. 


Appendix C. Influence of Bundle in High-Energy 
Muon Events 

High-energy muon events rarely consist of single par¬ 
ticles. Usually there is an accompanying bundle of low- 
energy muons, whose multiplicity depends on the pri¬ 
mary type and energy. It is possible to demonstrate that 
the influence of secondary particles on the leading muon 
energy reconstruction is negligible, and that information 
about the cosmic ray primary can be extracted using an 
additional observable. 

The accuracy of typical muon energy measurements 
can be increased by excluding exceptional catastrophic 
losses using the truncated mean of the energy deposition 
ll6^ . Since the high-energy muon energy estimate used 
in this paper relies only on the single strongest shower, 
the information used in the two reconstruction methods 
is fully independent. 

The approximate orthogonality of the two observ¬ 
ables can be demonstrated using only experimental data 
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by including information from the surface array IceTop. 
Since the leading muon rarely takes away more than 
10% of the primary cosmic ray energy, its presence has 
almost no influence on the surface size of the air shower. 
The signal registered by IceTop should therefore only be 
correlated with the properties of the cosmic ray primary. 

In Fig. C.l truncated mean and reconstructed muon 
surface energy are shown for the high-energy muon 
event sample as described in The lower two panels 
show the number of IceTop tanks registering a signal in 
coincidence with the air shower. The effect of varying 
the muon surface energy for a constant truncated mean 
is negligible, while in the inverse case a strong increase 
can be seen at the higher end. The result demonstrates 
that the total energy of the air shower, and consequently 
the size of the muon bundle, is not correlated with the 
measurement of the muon energy. On a qualitative level, 
it can also be seen that the truncated mean is related to 
the properties of the parent cosmic ray nucleus. 

For the quantitative interpretation of the truncated 
mean measurement, it is necessary to rely on simulated 
data, as illustrated in Fig. C.2 The true primary en¬ 
ergy distributions for proton and helium are clearly sep¬ 
arated. For the same nucleon energy, helium nuclei 
are four times more energetic than protons. The con¬ 
sequence is a substantially larger bundle multiplicity in 
the detector. To be distinguishable in the truncated mean 
observable, the energy deposition from the muon bun¬ 
dle needs to be comparable to that from leading muon. 
The relation between muon multiplicity and truncated 
mean is therefore less clear than in the muon multiplic¬ 
ity measurement as described in Section]^ 

A comparison between simulation and experimental 
data is shown in Fig. |C.3| The simulated curves are 
based on the simplified assumption of a straight power 
law primary spectrum. While a detailed analysis goes 
beyond the scope of this paper, the quantitative behav¬ 
ior of the experimental data conforms to the expectation 
that the average mass of the parent cosmic ray flux falls 
in between proton and helium. 


References 


[1] A. Karle [IceCube Collaboration], arXiv: 1401.4496 [astro- 
ph.HE]. 

[2] M. G. Aartsen et at. [IceCube Collaboration], Phys. Rev. Lett. 
Ill (2013) 021103 IarXiv: 1304.5356 [astro-ph.HE]]. 

[3] M. G. Aartsen et at. [IceCube Collaboration], Science 342 
(2013) 1242856 |arXiv:1311.5238 [astro-ph.HE]]. 

[4] M. G. Aartsen et ah [IceCube Collaboration], Phys. Rev. Lett. 
113 (2014) 101101 I arXiv: 1405.5303 [astro-ph.HE]]. 

[5] T. K. Gaisser, K. Jero, A. Karle and 1. van Santen, Phys. Rev. D 
90 (2014) 2, 023009 I arXiv: 1405.0525 [astro-ph.HE]]. 

[6] R. Abbasi et al. [IceCube Collaboration], Astropart. Phys. 42 
(2013) 15 I arXiv: 1207.3455 [astro-ph.HE]]. 

[7] R. Abbasi et al. [IceCube Collaboration], Astrophys. 1. 718 
(2010) L194 I arXiv: 1005.2960 [astro-ph.HE]]. 

[8] R. Abbasi et al. [IceCube Collaboration], Astrophys. J. 740 
(2011) 16 I arXiv: 1105.2326 [astro-ph.HE]]. 

[9] R. Abbasi et ah [IceCube Collaboration], Astrophys. J. 746 
(2012) 33 |arXiv:1109.1017 [hep-ex]]. 

[10] R. Abbasi et al. [IceCube Collaboration], Phys. Rev. D 87 
(2013) 012005 I arXiv: 1208.2979 [astro-ph.HE]]. 

[11] M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. D 89 
(2014) 102004 |arXiv:1305.6811 [astro-ph.HE]]. 

[12] https://web.ikp.kit.edu/corsika/ 

[13] T. K. Gaisser, T. Stanev and S. Tilav, Front. Phys. China 8 (2013) 
748 Iai-Xiv: 1303.3565 [astro-ph.HE]]. 

[14] T. K. Gaisser, Cosmic Ray and Particle Physics (Cambridge 
University Press, 1990). 

[15] M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. Lett. 
Ill (2013) 8, 081801 Iai-Xiv: 1305.3909 [hep-ex]]. 

[16] O. Adrian! et al. [PAMELA Collaboration], Science 332 (2011) 
69 I arXiv: 1103.4055 [astro-ph.HE]]. 

[17] M. Aguilar [AMS Collaboration], Phys. Rev. Lett. 114 (2015) 
17, 171103. 

[18] M. G. Aartsen et al. [IceCube and PINGU Collaborations], 
arXiv: 1306.5846 [astro-ph.IM], 

[19] H. S. Ahn, P. Allison, M. G. Bagliesi, J. J. Beatty, G. Bigongiari, 
J. T. Childers, N. B. Conklin and S. Coutu et al, Astrophys. J. 
714 (2010) L89 [arXiv: 1004.1123 [astro-ph.HE]]. 

[20] A. A. Kochanov, T. S. Sinegovskaya and S. I. Sinegovsky, As¬ 
tropart. Phys. 30 (2008) 219 |arXiv:0803.2943 [astro-ph]]. 

[21] T. Antoni etal. [KASCADE Collaboration], Astropart. Phys. 24 
(2005) 1 |astro-ph/0505413|. 

[22] J. Blumer, R. Engel and J. R. Horandel, Prog. Part. Nucl. Phys. 
63 (2009) 293 j arXiv:0904.0725 [astro-ph.HE]]. 

[23] D. R. Bergman and J. W. Belz, J. Phys. G 34 (2007) R359 
|arXiv:0704.3721 [astro-ph]]. 

[24] W. D. Apel et al. [Grande Collaboration], arXiv: 1206.3834 
[astro-ph.HE], 

[25] S. B. Shaulov, S. P. Beshchapov, K. V. Cherdyntseva, 
A. P. Chubenko, E. V. Danilova, Z. K. Zhanseitova, R. A. Nam 
and N. M. Nesterova et al, Nucl. Phys. Proc. Suppl. 196 (2009) 
187. 

[26] V. V. Prosin, S. F. Berezhnev, N. M. Budnev, A. Chiavassa, 
O. A. Chvalaev, O. A. Gress, A. N. Dyachok and S. N. Epi- 
makhov et al, Nucl. Instrum. Meth. A 756 (2014) 94. 

[27] M. G. Aartsen et al [IceCube Collaboration], Phys. Rev. D 88 
(2013) 4, 042004 jai-Xiv: 1307.3795 [astro-ph.HE]]. 

[28] W. D. Apel, J. C. Arteaga-Velzquez, K. Bekk, M. Bertaina, 
J. Blmer, H. Bozdog, I. M. Brancus and E. Cantoni et al, Phys. 
Rev. 0 87(2013)081101 |aiXiv:1304.7114 [astro-ph.HE]]. 

[29] W. D. Apel et al. [KASCADE-Grande Collaboration], Phys. 
Rev. Lett. 107 (2011) 171104 |arXiv:1107.5885 [astro-ph.HE]]. 

[30] W. D. Apel, J. C. Arteaga-Velzquez, K. Bekk, M. Bertaina, 


32 




J. Blmer, H. Bozdog, I. M. Brancus and E. Cantoni et ah, As- 
tropart. Phys. 47 (2013) 54 |arXiv: 1306.6283 [astro-ph.HE]]. 

[31] B. Peters 11 Nuovo Cimento 22(1961)4 

[32] S. V. Ter-Antonyan and L. S. Haroyan, hep-ex/0003006 

[33] J. R. Hoerandel, Astropart. Phys. 19 (2003) 193 |astro- 
ph/0210453|. 

[34] A. M. Hillas, J. Phys. G 31 (2005) R95. 

[35] V. 1. Zatsepin and N. V. Sokolskaya, Astron. Astrophys. 458 
(2006) 1 |astro-ph/0601475|. 

[36] V. Berezinsky, A. Z. Gazizov and S. 1. Grigorieva, Phys. Lett. B 
612 (2005) 147 |astro-ph/0502550|. 

[37] R. Aloisio, V. Berezinsky and A. Gazizov, Astropart. Phys. 34 
(2011) 620 |arXiv:0907.5194 [astro-ph.HE]]. 

[38] S. N. Boziev, JETP Lett. 54 (1991) 606 [Pisma Zh. Eksp. Teor. 
Fiz. 54 (1991) 603]. 

[39] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa and T. Sanuki, 
Phys. Rev. D 75 (2007) 043006 | astro-ph/06114181. 

[40] T. K. Gaisser, arXiv: 1303.1431 [hep-ph]. 

[41] A. Fedynitch, J. Becker Tjus and P. Desiati, Phys. Rev. D 86 
(2012) 114024 |arXiv:1206.6710 [astro-ph.HE]]. 

[42] P. Lipari, arXiv: 1308.2086 [astro-ph.HE]. 

[43] V. A. Bednyakov, M. A. Demichev, G. 1. Lykasov, T. Stavreva 
and M. Stockton, Phys. Lett. B 728 (2014) 602. 

[44] S. J. Brodsky, P. Hoyer, C. Peterson and N. Sakai, Phys. Lett. B 
93(1980) 451. 

[45] M. Britsch [LHCb Collaboration], Nucl. Phys. Proc. Suppl. 234 
(2013) 109. 

[46] E. Mountricha [ATLAS Collaboration], Nucl. Phys. Proc. Suppl. 
210-211(2011)37. 

[47] B. Abelev etal. [ALICE Collaboration], JHEP 1201 (2012) 128 
I arXiv: 1111.1553 [hep-ex]]. 

[48] B. Abelev etal. [ALICE Collaboration], JHEP 1207 (2012) 191 
I arXiv: 1205.4007 [hep-ex]]. 

[49] A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 97 
(2006) 252002 |hep-ex/06090101. 

[50] J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 94 
(2005)062301 |nucl-ex/0407006|. 

[51] M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason 
and G. Ridolfi, JHEP 1210 (2012) 137 |arXiv: 1205.6344 [hep- 
ph]]. 

[52] C. G. S. Costa, Astropart. Phys. 16 (2001) 193 |hep- 
ph/0010306|. 

[53] R. Enberg, M. H. Reno and I. Sarcevic, Phys. Rev. D 78 (2008) 
043005 |arXiv:0806.0418 [hep-ph]]. 

[54] J. 1. Illana, P. Lipari, M. Masip and D. Meloni, Astropart. Phys. 
34 (2011) 663 I arXiv: 1010.5084 [astro-ph.HE]]. 

[55] J. 1. Illana, M. Masip and D. Meloni, JCAP 0909 (2009) 008 
|arXiv:0907.1412 [hep-ph]]. 

[56] A. Fedynitch, R. Engel, T. K. Gaisser, F. Riehn and T. Stanev, 
arXiv: 1503.00544 [hep-ph]. 

[57] G. Gelmini, P. Gondolo and G. Varieschi, Phys. Rev. D 67 
(2003) 0173011 hep-ph/02091111. 

[58] S. 1. Sinegovsky, A. A. Kochanov, T. S. Sinegovskaya, A. Mis- 
aki and N. Takahashi, Int. J. Mod. Phys. A 25 (2010) 3733 
|arXiv:0906.3791 [astro-ph.HE]]. 

[59] M. Aglietta et al. [LVD Collaboration], Phys. Rev. D 60 (1999) 
112001 I hep-ex/99060211. 

[60] A. G. Bogdanov, R. P. Kokoulin, Y. F. Novoseltsev, 
R. V. Novoseltseva, V. B. Petkov and A. A. Petrukhin, Astropart. 
Phys. 36 (2012) 224 |arXiv:0911.1692 [astro-ph.HE]]. 

[61] M. G. Aartsen et al. [The IceCube Collaboration], Phys. Rev. D 
89 (2014) 062007 |arXiv:1311.7048 [astro-ph.HE]]. 

[62] M. G. Aartsen et al. [ IceCube Collaboration], arXiv: 1406.6757 
[astro-ph.HE]. 

[63] R. Abbasi et al. [IceCube Collaboration], Nucl. Instrum. Meth. 


A 703 (2013) 190 |arXiv: 1208.3430 [physics.data-an]]. 

[64] M. G. Aartsen et al. [IceCube Collaboration], JINST 9 (2014) 
P03009 |arXiv:1311.4767 [physics.ins-det]]. 

[65] E. I. Ahn, R. Engel, T. K. Gaisser, P. Lipari and T. Stanev, Phys. 
Rev. D 80 (2009) 094003 |arXiv:0906.4113 [hep-ph]]. 

[66] S. Ostapchenko, Phys. Rev. D 83 (2011) 014018 

|arXiv:1010.1869 [hep-ph]]. 

[67] K. Werner, T. Hirano, 1. Karpenko, T. Pierog, S. Porteboeuf, 
M. Bleicher and S. Haussler, Nucl. Phys. Proc. Suppl. 196 
(2009) 36. 

[68] J. H. Koehne, K. Frantzen, M. Schmitz, T. Fuchs, W. Rhode, 
D. Chirkin and J. Becker Tjus, Comput. Phys. Commun. 184 
(2013) 2070. 

[69] D. Chirkin and W. Rhode, hep-ph/0407075 

[70] D. Chirkin [IceCube Collaboration], Nucl. Instrum. Meth. A 725 
(2013) 141. 

[71] T. K. Gaisser, Astropart. Phys. 35 (2012) 801. 

[72] M. G. Aartsen et al. [IceCube Collaboration], Nucl. Instrum. 
Meth. A 711 (2013) 73 | arXiv:1301.5361 [astro-ph.IM]]. 

[73] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko and K. Werner, 
arXiv:1306.0121 [hep-ph], 

[74] M. Aglietta et al. [LVD Collaboration], Phys. Rev. D 58 (1998) 
092005 |hep-ex/980600I|. 

[75] J. Babson et al. [DUMAND Collaboration], Phys. Rev. D 42 
(1990) 3613. 

[76] V. N. Bakatanov, Y. F. Novoseltsev, R. V. Novoseltseva, 
A. M. Semenov and A. E. Chudakov, Sov. J. Nucl. Phys. 55 
(1992) 1169 [Yad. Fiz. 55 (1992) 2107]. 

[77] M. Ambrosio et al. [MACRO. Collaboration], Phys. Rev. D 52 
(1995)3793. 

[78] I. A. Belolaptikov et al. [BAIKAL Collaboration], Astropart. 
Phys. 7 (1997) 263. 

[79] E. Andres, P. Askebjer, S. W. Barwick, R. Bay, L. Bergstrom, 
A. Biron, J. Booth and A. Bouchta et al., Astropart. Phys. 13 
(2000) 1 |astro-ph/9906203 astro-ph/99062031. 

[80] G. Aggouras et al. [NESTOR Collaboration], Astropart. Phys. 
23 (2005) 377. 

[81] S. Aiello et al. [NEMO Collaboration], Astropart. Phys. 33 
(2010) 263 |arXiv:0910.1269 [astro-ph.IM]]. 

[82] J. A. Aguilar et al. [ ANTARES Collaboration], Astropart. Phys. 
34 (2010) 179 |arXiv:1007.1777 [astro-ph.HE]]. 

[83] P. Lipari, Astropart. Phys. 1 (1993) 399 |hep-ph/9307289|. 

[84] T. Adye, Proceedings of the PHYSTAT 2011 Workshop, CERN, 
Geneva, Switzerland, January 2011, CERN-2011-006, pp 313- 
318 I arXiv: 1105.1160 [physics.data-an]]. 

[85] A. Aab et al. [Pierre Auger Collaboration], arXiv: 1307.5059 
[astro-ph.HE]. 

[86] T. Abu-Zayyad et al. [Telescope Array and Pierre Auger Collab¬ 
orations], arXiv:1310.0647 [astro-ph.HE]. 

[87] P. Berghaus, T. Montaruli and J. Ranft, JCAP 0806 (2008) 003 
|ai-Xiv:0712.3089 [hep-ex]]. 

[88] M. G. Aartsen et al. [IceCube Collaboration], Astrophys. J. 809 
(2015) 1, 98 |arXiv:1507.03991 [astro-ph.HE]]. 

[89] M. V. Garzelli, S. Moch and G. Sigl, JHEP 1510 (2015) 115 
doi:10.1007/JHEP10(2015)115 |arXiv: 1507.01570 [hep-ph]]. 

[90] P. Desiati and T. K. Gaisser, Phys. Rev. Lett. 105 (2010) 121102 
I ai-Xiv: 1008.2211 [astro-ph.HE]]. 

[91] A. Margiotta [KM3NeT Collaboration], JINST 9 (2014) 
C04020 |arXiv:1408.1132 [astro-ph.IM]]. 

[92] R. Engel, A. Fedynitch, T. K. Gaisser, F. Riehn and T. Stanev, 
arXiv: 1502.06353 [hep-ph]. 

[93] P. Gondolo, G. Ingelman and M. Thunman, Astropart. Phys. 5 
(1996) 309 |hep-ph/9505417[ 


33 




Figure A.3: Top; Construction of differential energy 
deposition estimator. DOMs are represented by cir- p 

cles. The maximum lateral distance from the track up to | 

which individual data points are included in the recon- "g 

u 

struction can be varied depending on specific require- ^ 
ments. Bottom: Comparison between true and recon- 
structed energy loss in simulated event with parame- 
ters; iishower.reco = 1165 TeV (True Value; 852 TeV), 
cos 0zen,reco = 0.556 (True Value; 0.551) £^,1600 = 2493 
TeV (True Value; 1854 TeV). The shower energy corre¬ 
sponds to the highest single stochastic loss at approxi¬ 
mately 3000 m slant depth. Reconstructions using two 
different likelihood methods 1641 are shown for compar¬ 
ison. 



Figure A.4: Ratio between reconstructed and true 
shower energy for simulated events weighted to an ^ 
power-law primary cosmic ray flux spectrum. Around 
the peak the distribution can be closely approximated 
by a Gaussian distribution with a width varying between 
approximately 0.16 and 0.14. 



Figure B.l; Muon flux predictions from full shower 
CORSIKA simulation HTl and parametrization of the¬ 
oretical calculation ll53l . 
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Figure B.2: Effect of higher-order prompt flux correc¬ 
tion factor on all-sky muon flux derived from simulation 
using CORSIKA. The separation into cross section and 
re-interaction correction should be considered approxi¬ 
mate. 




Figure C.l: Top; Reconstructed muon surface energy 
and truncated mean ll6^ for experimental data. The 
sample corresponds to tracks with reconstructed angle 
within 37 degrees from zenith (cos 0zen > 0.8) in se¬ 
lection described in before exclusion of events with 
shower energies below 5 TeV. Red and blue boxes il¬ 
lustrate selection of data with approximately constant 
energy measurement. Middle: Number of IceTop tanks 
registering a signal in coincidence with muon track for 
fixed reconstructed muon surface energy (blue box). 
Bottom: Same for fixed truncated mean (red box). 
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Figure C.2: Parameter distributions separated by pri¬ 
mary cosmic ray type for simulated high-energy muon 
events with reconstructed surface energies between 30 
and 50 TeV. True primary energy (top) and muon bun¬ 
dle multiplicity at detector depth (bottom). 


Figure C.3: Top: Truncated Energy observable in COR- 
SIKA simulation weighted to GST-Global Fit flux and 
experimental data. Event selection criteria are the same 
as in Eig. |C.2| Bottom: Mean Truncated Energy ob¬ 
servable in dependence of reconstructed leading muon 
surface energy for simulated and experimental data. 
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